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ABSTRACT 


A dynamic model of 4 multicomponent-multistage distil- 
lation column with variable holdup and boundary lags is developed, 
and the programming for digital simulation is presented. The open 
loop transient responses for the mcdel undergoing feed composition 
disturbances and reflux changes are studied. These responses are 
then used for synthesizing transfer functions and closing the loop 
for digital control of the overhead product using a variable reflux 


rate. 
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Ll. Introduction. 


Consideration is given here to the dynamic performance of distil- 
lation columns with the ultimate object being column control studies 
by digital techniques. This is accomplished by mathematically model- 
ing a multicomponent-multistage column and using the model to generate 
column response. The column response is then used as an error source 
so that a control algorithm can generate a corrective action to move 
the system back to its desired control point. The model performance, 
error generation and control action are all done in a digital con- 
puter so that in essence direct digital control is being exercised. 

To have realistic control action, model responses should dupli- 
cate actual plant responses. In order to do this, it is necessary 
to include at the very least column hydraulics and holdup andidelay 
effects in the column end conditions. In this way, variable stage 
holdup is included in the model and the model will respond to changes 
in reflux flow which is one of the possible manipulable variables in 
this work. The work of Peiser and Grover [1] puts forth a model 
Similar to the one in this work. 

It is possible to generate model response without these secondary 
effects. Models of this type involve only mass transfer effects but 
can be quite complex if stage mixing and stage efficiency are incorpo- 
rated. In addition it is possible to have interaction between hydro- 
dynamics and mixing and efficiency effects. In the work presented here 
it is assumed that the stages are perfectly mixed and that the ef- 
ficiency is 100%. This in no way limits the model developed since 
allowance for both these effects could easily be incorporated into the 
model. 

Since a model will be of primary importance for response gener- 
ation, it is important to consider all possible physical processes 
which can effect a system response. Unfortunately, a model to do 
this becomes so mathematically unwieldy as to be useless. It has 
become the custom to make simplifying assumptions which produce a 
model that is relatively easy to work with. The sacrifice made is 
loss of realism in responses produced,iut this is not too important 


for purposes of study of general system behavior. However, to study 
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system control using models for response generation, it becomes 
necessary to put back in a model some of the "realism" lost by making 
the "usual simplifying assumptions." The problem of what to include 
in a model becomes a difficult one. 

Rosenbrock [2] in his excellent paper on distillation column 
control points out that the difference between observed and generated 
system response at small values of elapsed time is mainly due to 
secondary effects often not included in dynamic models. These second- 
ary effects are the hydrodynamic delays within a column and the delay 
and holdup effects in condenser and rebciler operation and external 
mass transportation. Unfortunately, it is at low values of elapsed 
time that errors for control action may be generated and so a model 
to be used for control studies should include these secondary effects. 

Since distillation has always been a much used operation in 
chemical engineering, the literature is voluminous regarding steady 
state analysis of multicomponent systems and dynamic analysis of 
binary systems. The literature on dynamic analysis of multicomponent 
Systems is somewhat less voluminous. In particular, work on multi- 
component systems wherein secondary effects are included is quite 
limited. Some recent publications in this area are the work of Sar- 
gent [3], that of Peiser and Grover [1], that of Moczek et al. [4] and 
finally that of Anisimov. [5] There are many other papers (see for 
instance [6], [7], [8]} besides these, but these are representative. 

As previously mentioned, then, the work presented here is the 
development cf a multistage-multicomporent model which allews for 
variable stage holdup, takes account of liquid hydraulics and, 
finally, allows for delays anc holdup in the condenser and reboiler. 
This model is used to generate responses to disturbances in the feed 
composition and the reflux flow. The responses are used in two ways 
as follows (1} to abstract the transfer function at any stage or the 


the error needed for control 


ct 
wD 


end of the column and (2) to genera 
action by means of a numerical algorithm. The transfer functions 
may be used for analog computer studies cf column control. 


One of the objectives of this work was to see if a numerical 


om 
(0 


don2 in reasonable time 


ia 


approach to column control stuiies could 


(9 


eted, the answer to this 


(0 


using digital computers. As might be exp 


is a function of the size of the system studied and the type of 
control involved, but indications are that work of this type can be 
done in reasonable time at reasonable cost. Another objective was 

to produce as accurate a model as possible with the criterion being 
duplication of actual plant response. Unfortunately, actual plant 
data was not available so that model checking could be done. Never- 
theless, the model as presented and used is accurate in the sense 

that at all times during a transient, the mole fractions sum to unity 
for relatively large values of the independent variable increment. 
This latter point is considered in some detail in the work of Mah 


et al. [9] and Sargent.({3] 
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2. Development of Model. 


General Model 





The above j a3 stage (j is the stage position index increasing 
vertically) shows all streams crossing the stage boundaries. All 
these streams with the exception of Q carry mass and energy. Q, the 
heat exchanged with the surroundings, is an energy stream only. 


A general dynamic mass balance for component "i" at stage "j" is 


a5 (We X ij BETS 5 Wy y, ); = (input - Youtput) ss (1) 


input = Loy) 4 + (LX, Dynlt cea + 7X J 


eo) eee Fae a 
[ SSF lLuooR SsF* Less 


houtput = lovy,), + ax,),| + [v Vi Ea eee, | 
iL) 2 SSP ivesp SSP Lissp 


A general dynamic energy balance for the same stage is 


ao“, h) EET 5 Wy HB), = (Yinput - pout) nerey (2) 
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Linput = low + on). )+ lF yey : Fhp, | + Vsor"yssr i Lge 


youtput = lov, + cua), | +Q+t+ LL + Laeghees| 


There will be (MeN + 2) equations of type (1) and (N + 2) equa- 
tions of type (2). The +2 accounts for the boundary condition equa- 
tions. In addition, there will be (MeN + 2) vapor-liquid equilibri- 


um relationships of the form 


/X (3) 


ee cies 


and (N + 2) mole fraction constraint equations 


Dy = ll X..=1 4 
3 oa5 7 a] (4) 


The equations as stated imply that (a) stage holdups, Wy and 
nae are perfectly mixed, (b) equilibrium between vapor and liquid is 
instantaneously realized at every part of a perfectly mixed stage, 
(c) holdup of vapor and liquid between stages is negligible, i.e. 
downcomer effects are not included, and (d) time to attain fluid 


equilibrium is small compared with that for mass transfer. 


Specific Model 

The general model, while reasonably complete, should be modified 
to reflect only the amount of "realism" required to accomplish the 
purpose of the model. Otherwise, the equations can become unwieldy 
and needlessly complex. With this in mind, the following additional 
assumptions are made about column operation: 

(a) the column operates adiabatically. (Q = 0) 

(b) vapor holdup is negligible compared to liquid holdup up to 
moderate pressures and can be neglected. 

(c) interaction between pressure drop changes, liquid hydraulics 
and mass transfer efficiencies can be neglected. 

These assumptions can be relaxed and the model modified to in- 
corporate any of these effects should it be desirable to do so. 

Equations (1) and (2) are next rendered dimensionless by dividing 


through by the feed rate, F. A dimensionless time unit, Tt, is 
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introduced, defined by the equation 
T= (FM, 8 (5) 


where the time, measured in units of tau, is equivalent to the number 
of tray turnovers based on the column feed rate. See Appendix VI for 
7+, 9 conversion factor. Further, since it is possible for ms to be 
different at every stage and vary with time, a particular Ws must be 
chosen for use with the definition equation (5). A reference stage 

is arbitrarily chosen and the holdup value on that stage at some point 
in time is designated WW ss and used in equation (5). Consistent 
with the above manipulations, it is also convenient to define a hold- 


up ratio 


Byac= “i i Ss C2) 


Thus, the model equations (1) and (2) now become 


a. ; 

aii .. al oss z : 

i, or a (XSinputs uoutputs) 65 See (7) 
dh 

gi eh, = b [coinput - Ran. | 
ae Ly (Linput OEE enerey nee (8) 


Variation in liquid holdup is allowed for by using expressions 
describing stage hydraulic behavior. Several detailed and very good 
correlations describing such behavior have been presented recently, 
i.e. see papers by Waggoner [16], Tetlow [17] and Peiser and Grover. 
[1] A modified version of the Peiser and Grover equation was adapted 


for use here as it provided an adequate but simple expression of the 





form ie = SOR See Appendix I for derivation. The equations used 
are 
ss, & 
.= 9p , + 30 eh 9) 
ag Pg * PPL GD 
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Ryy = yy + + fo, 9° 2h, a Pray 3°15 ©) (9a) 


An overall mass balance yields the following additional equation 


e 
Ry a de BCE DHE aes overall (10) 


Several additional relations are needed to account for the coef- 
ficients of equations (7) and (8) varying with time if any numeri- 
cal solution of the model equations is to be attempted. These are 
stated below. 


The equilibrium relations are assumed to be of the form 


Ai 


Kj = exp 7 B. + ia (11) 


Normally ae is a function of temperature, pressure and composition. 
In equation (11) it is assumed that there is no composition de- 
pendence and pressure drop through the column is negligible. 


The individual component enthalpies are given by 


hy eae, tab, (12) 


where the constants a. and b. are specific for each component. This 
same form of temperature dependence is also assumed to hold for 
individual component vapor enthalpies, Hs employing different 
constants. The values for these constants were obtained by cross- 
plotting data from Maxwell. [20] Consistent with the preceding 
relationships for enthalpies, the enthalpy of the liquid on stage 

j, designated nas is considered a function of composition and 
temperature, viz. 


h. = £(X; > T3) (13) 


j j 


For other than steady state conditions, Kay and rT are functions of 


17 


time so that the derivative of h, Can be written 


a 
dh, oe Skee) eh ed T 
meee hh = 5 |) ae (14) 
dQ jt (aX, 48 aT, do 


Assuming that ay can be expressed as a linear combination of 


individual component enthalpies, namely 


he = vh.X,, 15 
ee aes el 


it is possible to evaluate the partial derivatives in equation (14) 


and substitute back giving 


. Hu) fon 
et a cae pay, do 1S) 


Using equation (3) in the form 
= 17 
H(K,X i) 1 (17) 


and taking the time derivative gives 


dK. . 

. ee ee 18 
ehh SSS yl ods : 
where 
ae ei) pea) 
dona ij)| vg dg’? + &sSae (19) 


Substituting (19) and (18) into (16) and converting to dimensionless 


time, +, gives 
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° : pk; aX 
h, = Y(a,T, + b,)X + ja,X, (een (20) 
j ¢ ae | i 14 i 3 PK, X 5(Ay/T . cs 


The changes in composition from the integration of the mass 
balance equations (7) across At give values of Ki and a at 
(1 + Sr). From the X,., the T, at (1t + Ar) are obtained by a bubble 


point calculation. “nee equation (20) provides an independent 
estimate of fe since all of the terms are known at (7 + Ar). 
Equation (20) used in conjunction with equations (8), (9) and (10) 
suggests an iterative scheme to correct internal column flows ? 
across time. 

In summary, then, equations (4), (7) through (12), (15), and 
(20) are the model equations used for a numerical solution of the 
general stage dynamic response. 

It might be noted here that the individual component mass balance 
equations(7) can be considered to be a finite difference approxi- 
mation to (MeN + 2) coupled, nonlinear, partial differential equations 
in ee the independent variables being time and distance along the 
column. These partial differential equations are of second order in 
distance along the column and first order in time and also contain 


some first order terms in distance. This would indicate that the 


solutions do indeed describe traveling waves of composition changes. 
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3. Numerical Solution. 


Preliminary Discussion 

The data and results of this study were obtained using a CDC 
1604 digital computer to simulate the column. This section will treat 
the numerical integration of the general stage equations. The actual 
programs used for various boundary conditions are given in Appendix 
It. 

Before it is possible to simulate a transient response, it is 
necessary to know the initial steady state conditions for the nu- 
merical integration. All initial and final steady state conditions 
were calculated using Program I given by Hanson et al. [11] Results 
for the unperturbed and perturbed steady state conditions are tabu- 
lated in Appendix VII. Also, since the constants in the hydraulic 
equations (9) contain specific column parameters, some column design 
parameters for a typical column must be specified. See Appendix III 


for values used in this work. 


General Numerical Scheme 

A general approach for obtaining a numerical solution to a system 
of nonlinear differential equations with time varying coefficients is 
to set up a trial and error sequence for the integration with correct- 
ive iteration of the coefficients eventually forcing the equations to 
be satisfied. As was mentioned in the section on the specific model, 
equations (8) and (20) suggest an iterative scheme for correcting 
internal column flows within each time increment. Following this 
approach, the logic sequence used in programming a numerical solution 
for the response is as follows; 

(a) one stage is designated as the starting point for beginning 
the integration. 

(b) equations (7) for that stage are integrated across one Art 
using the Runge-Kutta-Gill algorithm as presented by Gill. [10] The 
coefficients are held constant at their last known values. This 
moves the ree to the point (7 + Art). 

(c) a bubble point calculation using equations (3), (4) and (11) 
is performed to establish a new T 


(d) liquid and vapor stream enthalpies are calculated from 
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equations (12) and (15). 
(e) equations (8) and (20) are evaluated. 


(f£) the test |n < ¢ (an arbitrary small 


j eqn(20) - i eqn(8 
number) is now used for convergence within a single Ar. If the test 
fails, for the first se citalat only, Ry is seme ce by an arbi- 
trary amount, AR 5? and Ray is evaluated as os aml L, is calcu- 
lated from equation (9), ue is calculated from equation (10), and 
steps (b) through (f) are repeated. 

(g) for the second and subsequent iterations, the last two known 
values of Raj are used to predict by the method of "reguli falsi" 
(21) that value of Raj which will cause the condition in (f) to be 
satisfied. When this condition is met, it is possible to proceed 
to the next stage and so on up and down the column until all stages 
have been integrated. 

(h) increment At and begin sequence again at (a). A flow diagram 


of this procedure is given in Appendix IV. 


Boundary Conditions 

In general, some control must be exercised or the column will 
not operate. For initial investigations it was assumed that the feed 
rate is ‘low controlled and the feed is saturated liquid at its 
bubble point; also, top product, bottom product and reflux stream 
flows were set and held. Thus, any changes in internal column flows 
are absorbed as holdup changes entirely. In a later section when 
control is discussed, these restrictions are modified. But, for 
the present, these assumptions are satisfactory for checking model 


accuracy and looking at the open loop response. 


lenis was found to be a good approximation after comparing results 
using equation (9a) and using an arbitrary ow . Use of equation (9a) 
requires a more complex iteration scheme whic Jis very sensitive to 
changes in internal flows and fails to converge except at very small 
increments of Ar. 


ee steps (f) and (g) several other similar schemes were tried 
which involved equating equation (8) to equation (20) and solving 
For... Ry and , were then calculated directly using the new 
valueJof£ ud All Sdch attempts failed to produce convergence within 
a single J Ar. 
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A total condenser was used in all simulation runs, both with 
and without an accumulator. Thus, reflux liquid is at its saturation 
temperature at all times. When an accumulator was incorporated, it 
was assumed to be an adiabatic mixer undergoing composition, temper- 
ature and holdup changes which in turn specify reflux stream con- 
ditions. 

The reboiler was treated as a large capacity equilibrium stage 
with a heat exchanger. No heat exchanger dynamics were included, 
but can be incorporated at any time. The reboiler duty was treated 
in two ways. 

(a) Reboiler duty (Q,) was set and held at the new known steady 
State value during the entire transient response. 

(b) Qe was considered a variable and allowed to float within a Art 
increment toward a new value which would correspond to the new equi- 
librium conditions in the reboiler at that point in time. 

The first situation was adopted for runs presented here. Physi- 
cally, this means that the controller set point for energy input to 


the reboiler is not changed. 


Delays 

A major use of the model response will be to derive transfer 
functions for control studies using the analog computer. In the 
many possible control schemes, it may well be that errors are detected 
at small values of elapsed time. Any time delay in the fluid trans- 
port at the column ends will affect the derived transfer function, so 
the model should have transport lag effects incorporated in it. One 
point where these lags occur are in the vapor flow from the top stage 
through the condenser-accumulator and reflux return line. Another is 
in the vapor flow in the reboiler return line. 

To introduce delays at these two points, a simple queing tech- 
nique was used. A plug flow situation was assumed in which the 
entering vapor stream variables Ws Tiss Hy ¥55 are those Leaving 
the top plate or reboiler at a particular increment in time and that 
these variables can be moved across the delay interval without mixing. 
Thus, the variables leaving the delay interval correspond to con- 


ditions 'n" Ar's ago. Diagrammatically, this can be represented as 
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follows: 
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4, Numerical Results 


Column Response to a Feed Perturbation 
To check the model, computer runs were made with a perturbed 
feed and variable end conditions. See Fig. 1 for physical system. 


The feed entering plate 4 consists of: 


Hydrocarbon Original Perturbed 
Composition Composition 

C4 ~22 +22 

C, 22 . 30 

Ce Aya .20 

Ce s23 225 


Variable boundary condition employed: 

(a) model with total condenser, constant Roa no accumulator 
and no fetayen : 

(b) model as in (a) with variable Raj added. 

(c) model as in (b) with accumulator added. 

_(d) model as in (c) with condenser and reboiler delays added. 

Moczek, Otto and Williams- [4] reported a substantial effect of 
the computation time increment used on the accuracy of the transient 
response; this was reflected to a large extent in their derived 
approximate transfer functions. RosenbEoe also discussed this 


point earlier. The main idea is the importance of using a compu- 


For a more complete discussion of this model, refer to a paper 
by Duffin, J. H., "Improving the Accuracy of the Dynamic Response 
of a Multistage-multicomponent Column Model," AIChE Symposium on 


systems and Process Control, 56th National Meeting, San Francisco, 
California, May, 1965. 


Z 
Moczek, J. S., Otto, R. E. and Williams, T. J., "Approximation 
Models for the Dynamic Response of Large Distillation Columns," 


AIChE Symposium on Process Control and Applied Mathematics, Vol. 61, 
1965, pp. 136-146. 


“Poscuntoce: H. H., Brit. Chem. Eng., Vol. 3, 1958, 
pp. 364-367, 432-435, 491-494. 
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FIG. 2 


8 STAGE COLUMN, TOTAL CONDENSER, 
EQUILIBRIUM REBOILER 
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4. Numerical Results 


Column Response to a Feed Perturbation 
To check the model, computer runs were made with a perturbed 
feed and variable end conditions. See Fig. 1 for physical system. 


The feed entering plate 4 consists of; 


Hydrocarbon Original Perturbed 
Composition Composition 
C, Ay ae | ihe 
C, «29 «30 
Ce ny ss) «20 
Ce Pe wen 


Variable boundary condition employed: 

(a) model with total condenser, constant Ray? no accumulator 
and no dois 

(b) model as in (a) with variable Raj added. 

(c) model as in (b) with accumulator added. 

_(d) model as in (c) with condenser and reboiler delays added. 

Moczek, Otto and Williams? [4] reported a substantial effect of 
the computation time increment used on the accuracy of the transient 
response; this was reflected to a large extent in their derived 
approximate transfer functions. Reseturoee” also discussed this 


point earlier. The main idea is the importance of using a compu- 


TOE a more complete discussion of this model, refer to a paper 
by Duffin, J. H., "Improving the Accuracy of the Dynamic Response 
of a Multistage-multicomponent Column Model," AIChE Symposium on 
Systems and Process Control, 56th National Meeting, San Francisco, 
California, May, 1965. 


2 
Moczek, J. S., Otto, R. E. and Williams, T. J., "Approximation 
Models for the Dynamic Response of Large Distillation Columns," 


AIChE Symposium on Process Control and Applied Mathematics, Vol. 61, 
1965, pp. 136-146. 


sR oseHbrOcK. H. Ho; Brit. Chem. Eng., Vol. 3, 1958, 
pp. 364-367, 432-435, 491-494, 
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FIG.4 
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tation interval of a certain minimum size even though computer time 


is correspondingly increased. The interval may be varied over the 
course of a computation to give desired accuracy and speed up solution 
time at larger values of elapsed time. 

For this study, the computation interval (41) was varied from 
.O1 to .3 (comparable to the .02 to 3.0 minute interval used by 
Moczek, Otto and Williams) and very little effect on accuracy was 
noticed. This was attributed to the small number of plates used 
in this model. See Fig. 2. In Fig. 2 the short, solid, vertical 
lines indicate the spread of points for Art varying between .03 and 
-30; the dotted lines join the data points for Ar = .1. The value 
of Ar = 0.3 was found to be the upper limit for a stable solution. 
A Ar of .1 was selected for use in all subsequent runs although 
computer time to attain steady state was on the order of 45 minutes with 
ample print out data for all ona No attempt was made to vary Art 
during a transient solution. 

Figures 3-8 are temperature, mass transfer, liquid flow and 
holdup response curves representing several stages in the column for 
varying column end conditions. These runs provided a check on the 
accuracy of the model in the serse that the model did move toward a known 
Steady state with the mole fraction sums always very close to unity 
and all time derivatives tended to zero within the error limit 
allowed by the iterative scheme. Table 1 lists representative values 
of these accuracy criteria over a typical solution interval. 

The effect of including tray hydraulics, lags and delays should 
be a Slowing of the mass transfer response with the final steady state 
being the same for all cases. Referring to Fig. 5, the temperature 
response of the top and bottom stages, and more specifically Fig. 3 
and 4,top stage and feed stage composition response for the C, 
fraction, it is evident that this is indeed the case. The responses 
flatten out and the approach to equilibrium becomes slower as the 


model is altered to include additional effects. An enlargement of 


S compater run time was substantially reduced by reprogramming 
in FORTRAN 63 vice 60. 
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the initial portion of Fig. 3 which is shown in Fig. 9 revealed some 
dead time except for the constant holdup model. With only eight 
Stages, dead time was very small and noticeable only three to four 
Stages away from the perturbation. For all intents and purposes, 
responses were immediate within the column. 

Referring to Fig. 6 and 7, the liquid flow responses on the top 
stage, feed stage and bottom stage were noticeably different from 
the mass transfer responses. The addition of hydraulics, an accumu- 
lator and then delays accelerated initial changes in the liquid flows. 
At larger values of elapsed Tt, the response curves crossed so that 
final approach to equilibrium was again slower as more effects were 
added to the model. It can be rationalized that, while the mass 
transfer responses show expected behavior, the lags introduced in the 
condenser-accumulator-reflux cycle and the speed of propagation of 
disturbance waves in liquid-vapor flows resulting from the inclusion 
of hydraulic interactions, cause changes in equilibrium conditions 
which badly upset the vapor-liquid flows initially. The rapid shift 
in internal flows should smooth out as the disturbance waves damp out. 
Also, the iterative scheme for movement of column variables implicit- 
ly alters vapor-liquid flows to meet the convergence test within each 
Ar. This produces rapid initial changes in flows which gradually 
smooth out as the time derivatives (hy, R, ,) pass through their 
maxima and become small. The final approach to equilibrium should be 
similar to the mass transfer responses. 

Representative effects of hydraulic equations on feed stage 
holdup are illustrated in Fig. 8. Table 2 gives the magnitude of 
holdup changes for all eight stages. It should be noted that care 
must be exercised in interpreting responses at very low elapsed t 
values. Any mismatch in the steady and unsteady state model will 
cause initial disturbances which must be allowed to die out before 


any desired disturbances is imposed on the system. 


Column Control Studies 
Control in the digital sense was done by letting the computer 
correct the responses through a control variable. This variable was 
adjusted via a control algorithm based on a generated error. The 
went Gay jeu L 1 Yeo Apes Niel LiL cbr cae woe tong “wie 
WET IAC VOCUEL Boyar ere 9 ‘mca Doc. sooner ees ( 
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control objective was to maintain constant C, 


composition in the 

overhead product by varying the reflux to compensate for disturbances. | 
Accordingly, the top product response to reflux changes was 

needed to determine the proper algebraic form of the control equation. 

Also, these open loop responses are necessary for synthesizing transfer 

functions for analog computer studies. Fig. 10 shows the actual model 

response to reflux perturbations. Fig. 11 is a plot of reflux rate 

versus overhead C, composition at final steady state. Results of 

these runs indicate that operation at a reflux rate around 1.0 limits 

the control action severely. In the vicinity of a 1.0 reflux rate, 

the steady state composition goes through a maximum. Thus, reflux 

changes would be ineffective in compensating for decreases in overhead 

composition at that point. As long as the perturbations into the 


system produce increases in overhead C, composition, then control can 


be maintained in this region. By Bee the initial operating point 
away from 1.0 reflux rate, fairly effective control can be realized. 

It was interesting to note also how fast the reflux disturbances 
were propagated through the column. Fig. 12 shows that there is a 
small delay before any change in liquid flow is detected in the bottom 
Stage and a slightly longer delay before any composition change occurs. 
This 4s expected. It takes a finite amount of time for disturbance 
waves to travel the column length. If the model is performing properly, 
with only eight stages the delay should be slight, but detectable. 

For this control study the model contained the following; 

(a) variable holdup 

(b) delays in condenser and reboiler lines 

(c) large capacity accumulator 

(d) liquid level control in accumulator and reboiler 
The computer was programmed to use a proportional plus integral 
control equation of the form r =r, + Ké + Ky Siedt 
Where the error, ¢, was generated by comparing the instantaneous 
value of C, in the accumulator to the initial starting value. An 
arbitrary error band of .001 in mole fraction was imposed, i.e. 
control action was taken if x - Kel > .001. X denotes composition 
of C, in the accumulator. Fig. 13 shows the results of varying K 


4 3 


and Ke The values Ke = K, = 1, gave an underdamped response with a 
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TABLE 2. | 
MOLE FRACTION SUM & _ DERIVATIVE VALUES ON 


FEED STAGE FOR A TYPICAL COMPUTER SOLUTION 
(At=0.1) 


MODEL WITH HYDRAULICS MODEL salt COLI AE & 
CUMULATOR 


vs tee *] tet] 
QNIS EQNIS 

90002 =| 189.9 9 -190.3 |i2xio"| | soeee |-100.9 -189. -199.9 | 

ee fl fw Ef 

ee i loc fo 
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; oe 

















TABLE ; 
INITIAL & FINAL VALUES OF HOLDUP RATIO 











REBOILER| | 


soo [oe [eon [ion ana nt | 
10.2100 .9822 | .997i | 1.0087 — .8913 |.9050 
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long settling time. Increasing the proportional gain to five, 
resulted in an overdamped response which settled out very quickly. 
However, a Kp = K_ = 10 produced an unstable response with increasing 
amplitude of oscillation. No attempts were made to optimize the 
response in any way. This portion of the work was merely to ascertain 
if control could be done in a reasonable time using digital control. 
Results were encouraging in this respect and the model performed 


adequately in all cases. 


ra} 


5. Transfer Function Synthesis. 


Technique Used 
The purpose of this section is to outline the technique used 
in synthesizing an approximate transfer function for the column and 
present the results of one specific case. 


The Laplace transform of X(t), 


on 


X(s) = Je Stxcepee (21) 
0 


can be transformed to the interval (0,1) by means of the substitution 


Then 
| 


X(s) = fet actos cjydr (22) 
e) 


This integral can be approximated by a Gaussian quadrature formula 


of order N appropriate to the interval (0,1), 


re s-1 
X(s) =X w,X(-log r,)r, (23) 
ist s* 2 


where (x,) are the roots of the shifted Legendre polynomials, 
# 


* 
Py (x) = Py Gt- 2x) 
; ; . 1 
and (wy) are the corresponding Christoffal weights. 


vTabulaced values of the (r,) and the (w,) are given for N up 
to 15 in Bellman, et al. [15]. *See also Beliman, et al. [12] for 
general approach to numerical evaluation of Laplace transforms. 
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Use of equation (23) in its present form restricts the X(t) 


values to a very small interval. The upper limit of the interval does 
not increase significantly as the order of N is increased, but the 
computational effort does. Therefore, to expand the X(t) interval by 
a factor, "a,'' one makes use of the multiplicative property of the 


transform and equation (23) becomes 


N 
B(siay, 2 % w,X(-a log r,)r et (24) 
a oe i’ i 


ae ae Be eer errr ere ae ere 


If the functional form of X(s) is known approximately and it is 
desired to improve the estimate of parameters contained in X(s), one 
can use a least squares curve fitting technique on the values of 
X(s) from equation (24) and the calculated values of the approiimate 


or assumed functional form of X(s), i.e. 


N 3 2 
S=2 Ace eqn(4) - X() assumed } 
Ss= 


where S is to be minimized with respect to the parameters in 


X(s) 


assumed" 

Curve fitting in the "s" domain in lieu of the time domain is 
preferred since the functional form of typical transfer functions 
for many complex processes involving heat and mass transfer have 
been investigated and reported im the literature. See for instance 
references (13], [18] and [19]. An alternative route for getting 
an approximate transfer function is to plot the normalized model 
response on semilog paper. The equation of the best asymptote to 
this curve yields the first term in the approximate equation of the 


-t/b -t/b 
t/b, +258 bibs ae a ras Next the 


response, i.e., X(t) = a,e 
difference between this asymptote and the actual response is plotted 
on semilog paper; the equation of the best asymptote gives the second 
term of the equation for X(t). In a like manner, the differences 


are plotted until as many terms as needed are obtained for the degree 
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of accuracy desired. A least squares fit could then be used as be- 


fore to optimize the parameters in X(s) with respect to 
approx. 
X(s) 


observed 


Transfer Function for Top Product Response to a Step in Feed 


The model contained an accumulator, delays in condenser and re- 
boiler, and liquid level controls in reboiler and accumulator. 


The transfer function G(s) for the situation 


STEP CHANGE C, (s) c (s) 
G 
IN FEED _—4\""_, TOP PRODUCT 


was the case studied. G(s) was assumed to be of the form 


-T Ss 
A eee Nene 
OE (eS Geer 
1 2 
and the values of K, Tay and T, were estimated graphically as 


described by Moczek et al. [13] from the open loop time responses. 


The results obtained were as follows: 


K v4 Ty) To 

1 ~3 tau 36.7 tau 4.1 tau 
or units units units 

1 -08 min 9.8 min 1.09 min 


If one denotes F(s) as the perturbation and X(s) as the response, 
the G(s) = X(s)/F(s). Thus, for a step in C, 
feed of 0.05, F({s) = .05/s and the functional form of X(s) to be 


composition of the 
used with equation (24) becomes 


-.38 
ee .05 e 
3(36.7s+1) (4.1841) 


A factor of 100 was used to expand the useful X(t) interval so 


that the actual form of X(s) used in the least squares fit was 
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X(s/100) 


100 SoS) Weta ws uae N 


The least squares program used to optimize the three parameters 


Oe ai? 


meter fit were K = .0536, Ties 36.2004 and Le 4.4585 or 


To) is given in Appendix V. The results for a three para- 


G(s) = —-b072 emacs 
S) ~ (36. 2004841) (4.4585s+1) 


(25) 
Using equation (25), X(s) = .05/sG(s) was inverted back to the time 
domain and plotted along with the model response for X(t). The 
results are shown in Fig. 14. The fit is excellent with only small 
deviations at large values of elapsed time--well beyond the control 


region. 
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6. Summary and Recommendations. 


A dynamic model has been developed which performs exceedingly 
well under the assumptions made in the mathematical development. 
Additional effects can easily be incorporated into the present model 
once the mathematics are specified. Preliminary control studies 
using the model for error generation indicate single loop and multi- 
loop studies can be done using a reasonable amount of computer time 
and appropriate transfer functions can be obtained from the model 
responses. All computation time can be substantially reduced by 
converting the present Fortran 60 programming to Fortran 63 or 
Fortran 4. Thus, the model can be a useful tool for design and con- 
trol work. 

Future work might be directed toward inclusion of pressure 
effects in the column, downcomer effects, reboiler dynamics and a 
partial condenser to extend the range of control variables and 
interactions which affect the system. Data on an actual operating 
column would be quite helpful in this respect, so that the model 
could be tailored to the actual plant conditions. In this way control 
studies and model responses could be correlated and checked against 
an actual operational plant in hopes of improving the design and 


control loops. 
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APPENDIX I 


DEVELOPMENT OF HYDRAULIC EQUATIONS 


The desired functional form of the expression relating liquid 
holdup to column flows is laa = £(Lj). The following development 
is adapted from an analysis for perforated plates done by Peiser and 


Grover [1]. 





weir height (ft) 

. = liquid height on the tray(ft) 

% | = liquid height in downcomer(ft) : 
AD; = cross sectional area of SORE SOS SNe ) 
A,,. = cross sectional area of tray(ft ) 

= liquid molar density(moles/ft>) 


h = height of the weir crest(ft) 


— 
tif 


m<— Mu. 
Ul 


width of the weir(ft) 


gravitational constant (£t/sce") 


mM! 


molar vapor flow(moles/time) 


molar liquid flow(moles/time) 


iM! 


+t 


c = correlation constant in the Francis weir equation 
Fp 5 = liquid flow resistance in downcomer 


Rj = liquid flow reversal resistance in downcomer 


Op = pressure drop for vapor flow and liquid head on tray abobe 
Up j = liquid velocity in downcomer 
oi = locally defined constant 


area under downcomer 


vapor velocity at tray (j + 1) 


perforation area 


total number of perforations at tray (j + 1) 


Hil 


orifice discharging coefficient (approx. 0.6) 


molecular weight of liquid 


molecular weight of vapor 
Using the Francis weir formula for large straight weirs 
B/z 
L, = kh (1) 
where 
k = c¢/Z7e 1. 
ee atn 


From the diagram 


h = Pry = Aaj (2) 


Using a pressure balance 


Hj 7 Pry = yt FRG t My (3) 
~ UD , 2 
eM) 0 cawer ar = Cr hel (4) 
Cy (2g) 


a2 


where 


Lid 2 


C..= ) 


- A 
Lj” 736(2g) *°L5 “v4 
In a similar manner 


2 


a vi Vi 4 
Vi, 
vats ang j+1 
sya ee DF” 
Therefore, 


re Bae 
Hoy vad eget 


where 
M, 2 
Cy jel ~ cz}, / (oR) 541 
Finally 
2 2 
Fg Py Fug yen © Sug s Fry 
1 Ws 


a | 4 
Pep Aaj Bie 14 


1 
where k “(k)¥ 
Liquid holdup, 


Lie = Cop Ag Pp) 2 (oD %) 5 


be 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


substituting equations (8) and (9) into equation (10) yields 
(11) 


| 1 45] 2 2 
= A ¢ cf) 72 a ¢ ‘ if V 
Weg 7 Pug gly * SPL PL Ee * Praha * Susy + Supe ] 


This expression is for a perforated plate. To simplify the 
analysis for a bubble cap tray, it will be assumed that 

(a) volumetric content of hydraulic gradient and downcomer 
equa{s volumetric content of risers and caps. 

(b) unidirectional flow predominates. 


(c) Francis weir equation holds. 


Let Pp j = effective liquid depth on plate j and A = effective 
surface area over which Pp § applies. Then 


we = oA ad (12) 
Substituting equations (1) and (2) into (12) yields 


fs 
Wes oy Alay + (L,/«) ] (13) 


Converting all flows to flow/unit feed flow and using a holdup 


ratio, Raj = weet dear equation (13) becomes 


Lr ss 
+ =. 

= 3 

Raj PF + 3p i (14) 
where 
2/s 
A A F 

veh a ite 

(Wiss (Wess & ly 


Differentiating equation (14) yields 


L u 4 
e 3 “3 Ld Je a e 
Ry = py + com Cs L,) +34 J? (on, | (15) 
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FORTRAN 60 PROGRAMS USED TO SIMULATE COLUMN RESPONSES 


APPENDIX II 


ON A 


CDC-1604 DIGITAL COMPUTER 


Nomenclature pertinent to all programs: 


X 


XIF 


Tif 


XISSF 


YISSF 


TEMP 
VAP 

QUID 

FL 

FV 

SSFL 

SSPV 

VAPMOH 

QUIMOH 

HFL, HFV, HSSFL, HSSFV 
RW 

RWDOT 


ENTHK, ENTHL 


ENTHV, ENTHW 


single component mole fraction in stage liquid 


single component mole fraction in entering 


liquid feed (Fi) 
fraction 


single component mole in entering 


vapor feed (F 
liquid side 


single component mole fraction in 


stream feed (Loop) 
fraction in 


component mole vapor side 


feed Voge) 


single 
stream 


temperature of stage liquid 

vapor flow 

liquid flow 

liquid feed entering a stage 

vapor feed entering a stage 

liquid side stream product leaving a stage 
vapor side stream product leaving a stage 
vapor molar enthalpy 

liquid molar enthalpy 

enthalpies of streams indicated 

liquid holdup on a stage 

liquid holdup derivative 


constants in the liquid enthalpy equation 
(h, = a,T, + b,) 
z 1] i 


constants in the vapor enthalpy equation 
(H, = a.T, + b,) 
1 4 i 
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A, B, C 


ALPHA, BETA, GAMMA 
DEL TAU 


BPERR 
DELERR 


RHOQD 


QC, QR 
NC, NS, JFEED, NCNTROL 


WEIR 


AREA 
FD 

CT 

GR 
RHOL 
HT 
WLRSS 
ALPH, BET 
DERIVM 
DERIVH 
QUNBAL 
NTALLY 


constants in the equilibrium equation 
(K, = exp(A,/T, Pee c,T,) 


constants in the Runge-Kutta-Gill integration 
computation interval for integration 


error limit allowed in bubble point calculation 


error limit allowed in convergence test 
([h ean (20) : H ,ean(8)I <e) 
individual component densities, Pra 
condenser and reboiler duties 


number of components, number of stages, feed 
Stage, number of increments over which inte- 
gration is to proceed 


weir length 

weir height 

effective tray area 

feed rate(F) 

constant, c, in the Francis weir equation 
acceleration of gravity, g2. 

total liquid density 

effective height of liquid on a tray 

WD ss 


constants q@ and g in equation (9) 
e@ 


ay 
Le from equation (20) 


i from equation (8) 


number of Ar intervals elapsed 


The following programs are for the model with varying boundary 


conditions. All programs utilize a total condenser and the hydraulic 


interactions, equations (9). 
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C THIS PROGRAM HAS AN ACCUMULATOR» BUT NO CON. 


PROGRAM DYNAMIC 


DIMENSION X(10950) 9TX(10950) »XIF(10950) sYIF (10950) sXISSFI10-50) YI 


1SSF (10950) sENTHK(10) s9EKC(10) » 


Fs ENTHL (10) sENTHV(10) sENTHW(10) 9EK (10) sALPHA(3) sBETA(3) sGAMMA( 3) » 
3DERIVM(10) » TEMP (50) »VAP(50) sQUID(50) »FL(50) o FV(50)sSSFL(50) eSSFV(5 
40) »SSPL(50) »SSPV(50) »SUM( 50) »VAPMOH(50) sQUIMOH(50) sDERIVH( 50) sQUNB 
5AL (50) sHFL (59) sHFV(50) sHSSFL(50) sHSSFV(50) #RW(50) sRWDOT(50) » 
6RHOQD( 10) sCON( 10) sFEED(10) sQUIDX(10) 2A(10)2B(10)9C(10) sDLYX(10) 


COMMON QUIDXsAsB9C 


OR REBOIL.» 


EQUILKF (A9BeCoT)tEXPF (A/(T+460.0)4+B4+C#(T+46020)) 


THALPF IY sZeT)2Y*T+Z 
READ 2eNCeoNSeJFEED »sNCNTROL 
2 FORMAT(4I15) 


READ 333s(ALPHA(1) sBETAC(I) sGAMMA(I) sI=193) 


333 FORMAT(6F12.6) 
1 FORMAT (8F 104) 
JLIM=NS+2 
READ lol (X(T oJ) sI=leoNC) sJ=lsJILIM) 
READ lol (XIF CI oJ) e¢IT=1lsNC) »J=29NS) 
READ Lol (YIFCI9J) sT=1l9NC) sJ=29NS) 
READ lel (XISSF(1sJ)sT2=1e9NC) sJe29NS) 
READ lel lYISSF( Il sJ)sT=leNC) sJ=2sNS) 
READ le(TEMP(J) sJ=leoJLIM) 
READ le (VAP(J)sJ=leJLIM) 
READ 1 e(QUID(J) sJ2=1lsJLIM) 
READ lLelFL(J) »J=2eNS) 
READ lelFV(JS) »J=29NS) 
READ 1o(SSFL(J) eJ=2eNS) 
READ leo(SSFV(J) sJ=2eNS) 
READ 1s(SSPL(J) sJ=29NS) 
READ lol SSPV(J) »sJ=29NS) 
READ 1»(VAPMOH( J) »J=1leJLIM) 
READ 1e(QUIMOH( J) eJ=lsJLIM) 
READ ls (HFL (J) sJ=2 NS) 
READ le (HFV(J)»J=29NS) 
READ lel HSSFL(J)sJ=29NS) 
READ 1s(HSSFV(J)»J=29NS) 
READ leo (RWIS) sJ=2e IL IM) 
READ1 +> (RWOOT (J) sJ=29NS) 
READ le (ENTHK( I) eo l=15NC) 
READ 1e(ENTHL(I)eF=19NC) 
READ lel(ENTHV( I)» I=1 NC) 
READ lelENTHW( I) oI =1 NC) 
READ 1elAl( I) eo I=19NC) 
READ 1e(B( 1) eT=leNC) 
READ l1selC(I)eI2#l5NC) 
READ 1,s0,ELTAUsBPERR»eDELERR 
READ 1e(RHOGD(I) »I=1sNC) 


C OR INITIAL = 136872026 


C OR UNDER PERTURBED CONDITIONS IS NOW SET AT 13247239 


uh 


DELAYS 


100 


1000 


FSN 


445 


444 


QR = 13247223948193 

WEITR=4,62 

HW=0225 

AREA=21 49 

FD=0227626 

CT=0042 

GR=3262 

RHOL=0.0 

DO 100 T=19NC 

RHOL=RHOL + X(I,JFEED) *RHOOD(] ) 
CONST=CT*SQRTF(2-*GR)#WEIR 
DENOM=CONST *¥RHOL 

HT=HW + (FD/DENOM) *#.6667* (QUID( JFEED) ) #26667 
WLRSS=RHOL *AREA#HT 
ALPH=AREA*#HW/WLRSS 

BET=AREA/WLRSS* (FD/CONST ) ** 26667 
LABEL =8HALPH BET 

PRINT 6sLABEL sALPHsSET 
FORMAT (A8 99F 1226) 

DO 9 J#3eNS 

RHOL £020 

DO 13 JT#leNC 

RHOL=RHOL + X(I9J)*RHOQD(I ) 
RW(J)=ALPH*#RHOL + BET*#(RHOL*#*, 3333) *(QUID( J) ** 6667) 
PRINT 1000¢(RW( J) »sJ=1loJLIM) 

FORMAT( 15H HOLDUP RATIOS ¢/912F10e6) 
PHI = O«l*DELTAU 

NTALLY=0 

NPRINT=0 

NPRINT = O 


445 SEQUENCE IS FOR STAGES BELOW FEEDSTAGEs 444SEQUENCE IS ABOVE 


JPATH=1 
JVARY=JFEED 
JLOW=JFEED 
JHIGHENS 

NUM=+] 

GO TO 444 
JLOW=3 
JHIGH=JFEED—1 
JVARY=JHIGH 
NUM=~-1 

DO 26 M=JLOWs JHIGH 
J=JVARY 
JVARY=JVARY+NUM 
LTIMES=1 


THIS SECTION THRU FSN 12 IS R~-K-G INTEGRATION FOR A GENe STAGE 


44 


DO 11 I=l»eNC 


TX(LoJ~-1)=X0 195-1) 
TX(TsJ) EXT.) 
11 TX¢(LoJ+1) =X 12541) 
DO 12 I=leNC 
QUE =0.0 
EQK=EQUILKF(A( 1) 5800) 9C( 1) »TEMP( J) ) 
TERM=VAP(J) ¥EGK+QUID(J) 
CON(T) =(VAP(J-1])*EQUILKF(CACT) . 301) os CCl) s TEMP (4-1))*TX01+J-1))+QUID 
L(J+i)*TX( 1,541) 
FEEDC I) SFLU J) #XIF CLs JD4FVI SD FY TEC Le J)4+SSFLEI IS FX SSF ULL oe J14SSFVEII) FY 
LISSE Cis J) 
OUTPUT=TX([sJ)*(SSPL(J)+SSPV(45) *E9K) 
CAYS(TONCT) -TXOT es JI FC TERM+RWDOT(L)) +FE EO! 1) -OUTPUT) *DELTAU/RWI J) 
DO 5 K=l1l93 
TX(TsJ)=T XL 95} +GAMMA(K) *{CAY-QUE) 
QUE =ALPHA(K) *®CAY+BETA(K) #QUF 
5 CAY=( CONC IT I)-TX(LoI)* (TERM+RWDOT(J)) +FEED( I) -OUTPUT) *DELTAUsRW( J) 
12 TX(LeJ)=TX( 195) +CAY/0-0-QUE/3.0 
31 SUM(J)=020 
DO 7 T=lseNC 
7 SUM(I)=SUM(J)+TX( I 9J) 
DO 8 T=lseNC 
TX¢TeJSJV=TXOC 125) /SUM( J) 
8 QUIDX(1)=TX(19J) 
CALL BUBPT(TEMP(J) »NC»sBPERR) 
VAPMOH( J) =C.0 
QOUIMOH(J)=0-0 
DO 15 T=leNC 
EK (CT)=EQUILKF(A( 1) sBU1) os CCl)» TEMP(J)) 
DERIVM(T)=( CONC IT) -TX( 1 ed) FCEK( IL) *¥VAP(JS)4QUID(J)+RWDOT(J)) +FEED( I)- 
1 TXC1sJ)*(SSPL(I)4SSPV(J)*EK OCI) )/RWOCJ) 
QUIMOH( J) =QUIMOH(J)+THALPF (ENTHK( I) sENTHLC I) » TEMP( J) )*® TX 19d) 
15 VAPMOH( J) =VAPMOH(J)+THALPF(ENTHV( I) sENTHW( I)» TEMP( I) DFTXC Ts J) #FEK CI 
1) 
TERM1=0.20 
TERM2=020 
TERM3=0.20 
TERM4=020 
DO 21 IJ=1l»NC 
TERMIL=TERMI+(ENTHK (I) * TEMP (J)4FNTHL(OIT)) FPR ERIVMCOT) 
TERM2=TERM2+EK(1)#DERIVM(1) 
TERM3=TERM3+EK(C 1) *® TX 195) RCAC T)/C0C TEMP (I) 446020) #*¥2)-Cl(1)) 
21 TERM4=TERM4+ENTHK( IT) * TXT 95) 
DERIVH( J) =(TERM1+(TERM2*®TERM4) /TERM3) 
FD THALP =VAP ( J-1) ® VAPMOH(J-1)+QUID(J+1) *QUIMOH(J+]1)+FL( J) *HFL(J)+FV 
1(l) *HFV(J)+SSFL (JS) FHSSFLII)+SSFV(5) FHSSFV( J) 
GO TO (33933%221)sJPATH 
33 QUNBAL (J) =( FDTHALP-QUIMOH(J)*(SSPL(J)+RWDOT(J)+QUID( J) )-VAPMOH( J) * 
1(SSPV(J)4VAP(J))) /RWOJ) 


C THE NEXT STATEMENT COMPUTES THE DERIVATIVE DIFF FOR CONVERGe TEST 


DIFF=DERIVH(J)-QUNBAL(J) 
TF CABSF (DIFF )-DELERR) 25925922 


CC FSN 


22 


101 


2a 


35 


24 


a7 


25 
26 


Cc FomN 


132 


502 


504 


22 THRU 26 IS ADJUSTMENT OF COEFFe BY REGULI FALSI 


LTIMES=LTIMES-1 

RHOL=0e20 

DO 101 IT=1lsNC 

RHOL=RHOL+TX(19J)*RHOQD(1) 

EXP=RHOL**, 3333 

IF(LTIMES) 24923523 

RWO=RW( J) 

RWREF=RW( J) 

RW(J) = RW(J) + PHI 

RWOOT(J)=el 

QUID( J) =((RW(JS)-ALPH*#RHOL ) /(BET*EXP) )##1-5 
SUMFDS=VAP( J-1)+QUID(J41)4FL(5)+FV(JS)+SSFL(J)+SSFV(J) 
VAP (J) =SUMFDS-QUID(J)-RWDOT (J) 

DIFFO=DIFF 

GO TO 44 

SLOPE=(DIFF-OIFFO) /(RwW(J)—-RWO) 
RWNEW=RW(J)-DIFF/SLOPE 

RWDOT (J) =(RWNEW-RWREF ) /DELTAU 

QUID( J) =€ (RWNEW-ALPH*#RHOL )/(BET#EXP) ) #165 
VAP (J)=SUMFDS-QUID(J)-RWOOT(J) 

DIFFO=DIFF 

RWO=Rw({ J) 

RwW(J)=RWNEW 

IF (LTIMES+30) 37944944 

LABEL=8HEND TEST 

PRINT 6sLABEL 29 JsQUNBAL (J) sDERIVH( J) sQUID( J) s VAP (J) sRW( J) sRWDOT( J) 
GO TO 30 

DO 26 I=1sNC 

X(IsJ)=TX( 1 9J) 

JPATH=JPATH+] 

GO T0(199%9199%9201) sJPATH 


199 THRU 219 INTEGRATES THE CONDENSER/ACCUM SECTION 


J=ENS 

DLYV=VAP(J) 

DLYHV=VAPMOH (J) 

DLYHL=00 

DO 502 I=1sNC 

EK( IT P=EQUILKF(CA(T) sB( 1) 9CC1) sTEMP(J)) 
DLYX(I)=EK(C I) *XC19J) 

QUIDX(1)=OLYX(I) 

TEMPCON=TEMP (J) 

CALL BUBPT( TEMPCON »NC sBPERR) 

DO 504 [=1sNC 
DLYHL=DLYHL+THALPF (ENTHK( I) s9ENTHL( I) sTEMPCON) ®DLYX(T) 
CONDMOH=DLYHL 

Qc=DLYV #*(DLYHV -DLYHL ) 
CONLIQ=QUID(J+1)+QUID(J+2) 


LTIMES=] 
RWREF=RW(J+1) 
RWDOT (J+1)=DLYV -CONLIQ 
RW(J+1) =RWREF+RWDOT( J+] )*DEL TAU 
205 DO 206 I=1lsNC 
206 TX( IT oJ+l) =X 1 eI5+1) 
DO 207 I=1leNC 
QUE=0.0 
CAY=(DLYV#DLYX(1)-—(CONL 1Q+RWOOT (J+1))* TXT eJ+1)) #DELTAU/RW( J+1) 
DO 208 K=193 
TXC i sJt1LIPHTXC3 9 J4+1 )+GAMMA(K) *( CAY-QUE) 
QUE=ALPHA (K) *#CAY+BETA(K) #QUE 
208 CAY=(DLYV*DLYX(1)-—(CONLIQ+RWDOT(J+1))*TX(194+1))* DEL TAU/RW( J+1) 
207 TXCTsJ+1)=TXC 19 J+1)4+CAY/6e0 -— QUE/3.0 
SUM (J+1)=020 
DO 209 I=1¢sNC 
209 SUM(JS+1)=SUM0J+1)4TXO19d+1) 
DO 210 I=1eNC 
TX( 1 eJt+1l) =TX(19I4+1)/SUM( 541) 
210 QUIDX(I)=TX(1-J+1) 
CALL BUBPT(TEMP(J+1)} »sNCsBPERR) 
QUIMOH( J+1)=040 
DO 211 I=1lsNC 
EKC(I)SEQUILKFCACT) so BC I) »CCI Ds TEMP (J+1) ) 
DERIVM( I) =(DLYV#DLYX(1)-(CONL IQ+RWDOT ( J+1) VF TXC Ts J+1))/RWOJS+1) 
211 QUIMOH(J+1) =QUIMOH(J+1)4THALPF (CENTHK (1) s9ENTHL (1) TEMP(J4+1))*TX0T oJ 
lea 
TERM1=0-20 
TERM2=0-0 
TERM3=020 
TERM4=0,0 
DO 212 I=lsNC 
TERMI=TERM] +(ENTHK( I) * TEMP ( J+] )+ENTHL (1) )*#DERIVM( 1) 
TERM2=TERM2 +EKC(1)*DERIVM(1) 
TERM32TERM34EKC( 1) FTXOCToStLIFCACT) SCC TEMP( S41) +4600) ##2)-Cl1)) 
212 TERMS=TERM4+ENTHK( 1) *TXC1sJ5+1) 
DERIVH( J+1)=( TERM] +( TERM2*TERM4) /TERM3) 
QUNBAL (J+1)=(DLYV * CONDMOH- ( CONL 1Q4+RWDOT (J+1) ) *QUIMOH( J+1))/RWOJ 
I+1) 
216 DO 219 I=19NC 
219 XCTeoJ+1L)=TXC1sJ+1) 
GO TO 445 


FSN 201 THRU 228 INTEGRATES THE REBOILER SECTION 


201 J=2 
LTIMES=1] 
GO TO 44 
221 QUNBAL (J) =( FOTHALP-VAP (J) *VAPMOH(J)-QUIMOH(J)*(QUID(J)+RWDOT(J))+Q 
IR)/RWOJ) 
DIFF=DERIVH(J)-QUNBAL (J) 
IF (ABSF (DIFF )-DELERR) 22592259222 
222 LTIMES=LTIMES-1 
IF (L TIMES) 22492239223 


223 RwO=Rwi( J) 
RWREF=RW( J) 
RW(J)=RW( J) 40601 *DEL TAU 
RWDOT(J)=0-0]1 
VAP (J) =QUID(J+1)-QUID(J)-RWDOT (J) 
DIFFO=DIFF 
GO TO 44 
224 SLOPE=(DIFF-DIFFO)/(RW(J)-RWO) 
RWNEW=RW(J)-DIFF/SLOPE 
RWDOT (J) =(RWNEW-RWREF) /DELTAU 
VAP (J) =QUID(J+1)-QUID(J)-RWDOT( J) 
DIFFO=DIFF 
RwO=RW( J) 
Rw(J)=RWNEW 
IF (LTIMES4+30) 226944544 
226 PRINT 227 sQUNBAL(J)sDERIVH(J) 
227 FORMAT(13HREBOILER LOOP >2F208) 
GO. FO) 30 
225 DO 228 IT=1lsNC 
X(T sJ-1)=TXC19J) 
228 X(1leJ)=TX(1 oJ) 


C FSN253 THRU STOP 1S THE DATA PRINTOUT SECTION 


253 NTALLY=NTALLY+1 
IF (NTALLY-NPRINT)4917917 
17 PRINT l18eNTALLY »sDELTAU 
18 FORMAT(CI7HILINCREMENT NOe = I[5/34HOINDEPENDENT VARIABLE INCREMENT = 
1 F104) 
PRINT 19 
19 FORMAT (116HOSTG MOL FRA SUM H UNBALANCE H DERIVATIVE TEMPER 
LATURE LIQUID FLOW LIQ ENTHALPY VAPOR FLOW VAP ENTHALPY ) 
JSTG=0 
JL IM=NS+1 
DO 27 J=2sJLIM 
PRINT 209JSTGsSUM(J) sQUNBAL (J) ,DERIVH( J) » TEMP( J) sQUID( J) sQUIMOH(J) 
1 sVAP (J) sVAPMOH (J) 
2C FORMAT(1498E14-8) 
27 JSTG=JSTG+] 
PRINT 271+9CsQR 
271 FORMAT(18HOCONDENSER DUTY = E15¢e¢8/17HOREBOILER DUTY = E1508) 
JSTG=0 
DO 29 J=2sJLIM 
PRINT 28sJSTGs(X(1eJ) sT=leNC) »sQW( J) sRWDOT( J) 
28 FORMAT (29HO STAGE NOe = 13/(10F13410)) 
29 JSTG=JST1G6G+ 1 
NPRINT = NPRINT + 1 
IF (NTALLY-NCNTROL) 4930530 
30 STOP 77777 
END 
c SUBROUTINE FOR BUBBLE POINTS 


SUBROUTINE BUBPT(T sl »BPERR) 
DIMENSION GUIDX(10)9A(10) 9B(10) 9C(10) 


ae 


COMMON QUIDX sAsBoC 

EQUILKF (AsBeCesT)SEXPF(A/(1T+460.0)+B+C#(T4+46060)) 
KTIMES=1 

SUMY=0.20 

DO 4 T=leL 

Y=EQUILKF (ACT) sBC IT) eo CCl) eT) #QUIDX(T) 
SUMY=SUMY +Y 

IF (ABSF (SUMY-1.20)-BPERR)898s5 
KTIMES=KTIMES-1 

IF(KTIMES) 79696 

SUMY O= SUMY 

TO=T 

T=T+1020 

GO TO 3 
SLOPE=(SUMY-SUMYO)/(T-TO) 
TN=((120-SUMY ) /SLOPE )4+T 
SUMYO= SUMY 

TO=T 

T=TN 

GO TO 3 

RETURN 

END 

END 


C 


Pr Ce aot? VES 
PRUGRAY SY NAM Le 


Tals PROG n a4 MAS NO DEukeS Oe ACCUCULA TOR 


333 


WIGENS LON ACG0 55075 1420550), X2F (20S 0S VIF C107 50) 541556 Clos oa 
SSP C105 oO) EN TK Oo 

2 ENTHUCLO),ENTHVILO)ENTHW(10),2K(10),ALPRA(3),BETA(S),GAMMA(3), 
SDERIVHC10) + TEMP (50) .VEP (90), 0UID(56),FL(S0).FV C90), SSrLCe0 yp soem 
40).+SSPL(50)>SSPV(5U),SUM(50),VAPMOH(50),QUIMOH(50),DERIVH(50),QUNB 
PAL (30), 4FL(50),HFV(50),HSSFL(50),HSSFV(90),RW(50),RWD07T(50). 
6Re000(610),CON(10),FEED(10),@UIDA(10),A(10),B8(10),0(10) 
COMMON QUIDA,4,58,C 

CQUILAF (A, B.C, T)=EXPFCA/(T+460.0)+B+C=(T+460.0)) 

THREP TO ta 2 el im Ve eZ 

READ 2@,NC,NS;, JFEED»NCATROL 

FORMAT (415) 

READ S55, (ALPHA(1),BETA(I),GAMMA(1),121,53) 

PORM AOR tLe 20) 

PURMAd VOC 10440 

sLIM=NS<-2 

READ SOCK BJ Uy oe dG Vee Ti) 

READ Toke Ole l= Noe 2 ano 

RSA re a Cty ee Cig dd5 1 =i NC Ja25N5) 

READ 2; CCALSSr (la Jal 1TeNC) ye =25NS2 

READ Se CCT ISS (15575 |=25 NC) eJS2eN5) 

READ a7 CLEMO) deco LEM 

READ 29 (VAP CU) a JS1s JL IM) 

ReAD eee re I GIN a ars | ares melee a 

READ 2, (FL(J),J=2,NS) 

READ PEN Cy ae NS 

READ TVA SSE UDG e2NS) 

READ Le (SSEVCJ).Js25N5) 

READ WresSbewmd) pda2,.NS) 

READ ee (Oo Se VC) Wle es Noo 

READ 2+ ( VAPMORC I), gS o eck IM) 

KeAD 2. COUINOHC Jig Jet eg) 

RERD 2, (AFL( J), Je2sNS) 

Heed ves Car VC le2a NSD 

REL Se thost lt dose oho 

EAD 2, CHSSPV(CJ)> Je2. NS) 

READ 2s (RWCJ)s J=2,J5L18) 

READL,. (RXDOT(J),4=2,NS) 

ACAD gk Stee Le eGo 

READ Le (ENIHL(OI),1=1,NC) 

READ 21,¢ENTHV(OI),]1=1,NC) 

xCAD 1,¢(ENTHW(]),]=1,NC) 

READ ciara er “ae Ge Jars Uae fes tara 

READ 2,-(BC1),1=2,NC) 

22h) aN Otol ete NGO 

READ 1,D0EL.TAU,GPERR, DELERR 

EAD Lo (RWOOD¢ 1 yl ata NG) 


QC = 13247. 23948/93 


fei net. o2 

maeG.25 

AnRtA=22.49 

rD=0.27626 4 


100 


445 


444 


44 


a3 


tz 
od 


DG 106 j{=1,NC 

RHOL=RHOL + X(€1,JUFEED) *RHOOD( 1) 
CONST=CT«SGRTF (2. eGR) eWEIR 
DENGM=CONST¥RHOL 

HT=HW + (FD/DENOM) **.6667~(QUID(UFEED) )**.6667 
WLRSS=RHOL*AREA*HT 
ALPH=AREAWHW/ALRSS 
BET=AREA/WLERSSe(FD/CONST) we, 6667 
LABEL =8HALPH SET 

PRINT 6,LASEL,ALPH,BET 
FORMAT(AB,9F12.6) 

DO 9 J=3.,NS 

KHOL =G¢.0 

BO. 1s Le taane 

RHOL=RHOL + X(J,J)*RHCQD(]) 
RWCJSALPH*RHOL + BET#(RHOL**.3533) *#(QUID(J)**.6667) 
PRINT 2000.(RW(J),J=1,JLIM) 
FORMAT(15H HOLDUP RATIOS ,/,12F10.6) 
PHI = 0.1¥*DELTAU 

NTALLY=9 

NPRINT=0 

JPATHE=L 

JVARY=SJFEED 

JLOW=JFEED 

JHIGHENS 

NUMs+4 

GO TO 444 

JLOW=S 

JRIGH=JFEED-1 

JVARY=JHIGH 

NUM=-2 

00 26 M=JLOW, JHIGH, 

J=ZIVARY 

JVARY=JVARY*+NUM 

LTIMES=1 

DO 12 1=1,NC 

TXC1.J-7) =X 1,J=1) 

le 5 ee Ow & Oh aeew 

TXC1,J3¢1)5X(1,J541) 

DO 12 [21,NC 

QUE=0.0 
EQKSEQUILKF(CA(T),B(1),C(1),TEMP(J)) 
TERM=VA9(45) *EOK+QUID( LU) 

CONC ITI S(CVAP( J-1) eEQUILAF (ACI ),B°¢1),C001),TEMP(J-1))#TXC1,4U-1))*QUID 


LCJd*l) +i xClsJde1) 


FEEDCIDSFLAJ0 eX IP CIs J) eP VC) eV IP C1] ot +SSFUCU)eXISSE Cla Jd *+SSFVCU) eY 


LISS (i590 


OUTPUT=TX(1,J)2%(SSPL(J)+SSPV(J) *EOK) 
CAY=(CONCTI-“TX(I,U)*(TERM*+RWDOT(J)) +FEEDC I) -OUTPUT) «DELTAU/RW( JU) 
CI 5 K=1,.3 

TACT, SJVSTXCL, JS) + GAMMA (K) «(CAY-QUE) 

CUESALPHA(K) #CAY*357~ ACK) ©QUE 

CAY=CCONC TIM TXCT, gas TERM*RADOTCU)) +FEEDCI)-OUTPUT) *DELTAU/RW(Y) 
SACI, S02 S7X01,J5)+CAY/620-GUE7S.0 

SUM(J)=0.0 

vo 7 i=1.,NC 

SUM( J) =SUMC J) *TX( 1,5) 

CO 8 I=:.NC 

The ToS TKC 1s Jd 7SUNCI) 55 


& 


ale) 


igh 


33 


Ze 


1014 


23 


24 


QUIDACT)=HTXC] 
GALL SUBPT(TE 
VAP MOR C Joe On) 
SUMO GJS Ol cn 
DO 15 L=1,NC 

ex(lPSEOULERE CACTI) oot idwct ld. Veneta) 
DERIVMCTISCCONCIT)-TXC1L,45) ® (CEKC 1) we VAP(J)*+QUID(J)*RWDOT(U))+FEED(I)- 
TACi» J) *(SSPLOJ)+SSPV( J) wEK(1)))/RWOJ) 

QUIMOn (J) SGUIMOH(U)+TRALPF CENTAK(I),ENTHL(I),TEMP( J) )*®TX(1 0) 
VAPMOH(J)=VAPMOH( J) *TRALPFCENTAHV( I), ENTHW(T),TEMP(J))*TXCI,J) @EKCI 


J 
MP(J),NC,BPERR) 


TERMG=0.C 

TERMZSTERMS + (CENTHK (I) *® TEMP CU) +ENTHL OI) ) *DERIVM(I) 
TERM2=TERM2*+EK( I) *«DERIVMC(T) 
TERMS=TERMS*EKC 1) eTX(C1,J)e(ACT)/ (0 (TEMP (J) +460.0)**2)-C(1)) 
TERM4=TERM4*ENTHK( 1) *1TX(1,0) 

DERIVACJ)=(TERML+(TERM2*TERM4)/TERMS) 

FDTHALPSVAP( J-1) *VAPMOH(J-1)*QUID(J+1) *QUIMOH( J*4)*FL (J) #HFL (J) +FV 


1(J)*HFV(J)+SSFL( J) *HSSFL (J) +SSFV( J) *HSSFV( J) 


GO TO (33,335,221), JPATH 
QUNBAL( J) = (FDTHALP-QUIMOH( J) *(SSPL (J) *RWDOT( J) *QUID(C J) )-VAPMOH(J)* 


ICSsevCIy eV art Joy) /RW A) 


DIFF=DERIVH(U)-QUNBAL( J) 

IF CABSF (DIFF )-DELERR) 25,25, 22 
LTIMES=LTIMES-1 

RAQL=0.0 

DOF 1g PaiG.NC 

RHOLS=RHOL+TX(1,J)*RHOGD( I) 

EXP=RHOL**.33585 

TECLUIMES). 24,5233 28 

RwO=RW( J) 

RWREFZRW( J) 

Rv( J) = RWCJ) + PHI 

RWDOT(U)=.1 

QUID( J) =( (RWC J) -ALPH*RHOL)/(BET*EXP) ) "1.5 
SUMF DS=VAP(J-1)+QUID(U+1) 4+FL (JU) + F VC) +SSPFL (CU) +SSFV(J) 
VAP (J) =SUMFDS-QUID(U)-RWDOT( UY) 

DiFrO=DIFF 

SO TO 44 

SLOPE=A(DIFF-DIFFO)/(RW( J) -RWO) 
RwNEweRwW(J)-DIFF/SLOPE 

RWDOT( JU) = (RWNEW-RWREF )/DELTAU 

QUID( J) = ( (RWNEW-ALPH*RHOL)/(BET*EXP) )*#1,5 
VAP (J) =SUMFDS-QUID(J)-RWDOT (CU) 

DIFFO=DIFF 

RAO=RW( J) 

Rw J) SRwNEW 

IF(LTIMES+30) 37,44,44 

LABELZBHEND TEST 

PRINT 6,LABEL,J,OUNBAL(J),DERIVAH( JU), QUID( JU), VAP(U),RW(J),RWDOT( J) 
GG 0 36 

DO 267 1=1, NC 

a Get rgee at ie Ga 

JPATHESUJPATHTAL 

BO 10t 199,200,201 opera 

SSS 

QUID(J*2) =VAP(J)-QUID(J+1) 

DO 502 [=4,NC 


ee Pee API ae 2 ope dh 


>. THe 
i ci Ge 
4 ott 
Bes 
tH 
2 


Cy ee a 
a Ces 


Cree. 
-tu 


ee 


m 


“MPC deaisACesPenn) 
0 


1 T- 

oe 

-- C ent? 
CL. mse OC.) Us ee we 4) 


fy 
4m 4 Uw  € 


Acts ts 


(ENB a SER Go 
2< Ul 
oc © 
= & 

ws 
aoc 


as 
GC 


TMOnC Sele PaAL Pe CENT AA (Ci) ENTAL CI) TEMP C41) XC] e 


P( J) #CVAPMGAN( J) -CUIMSAC Je i) 


u 
“<= 


oA tw aogQgw~e COC Fr Ges wv 
e~ {3 oe 


36) '— C..G31>iF* 


en 


= 


ACP HVA (Sw VARPMOA( IU) -QUIMOH(U) © (QUID( J) *RWDOT( J) )+0 
tas ae 
ee Oe 
222 LTIMESSLTIMES@4 
Ir (LTimes)224, 220, 223 
225 RAOQ=zRW CU) 
PrREFFRHC J) 
Rvi( JIERWCJI*0.02*DELTAU 
RA#DITCII=5-072 
VAP (J) =QUI0(U+1)-QUID( J) -RWDOT (J) 
DIF FOSDIFF 
GO Td 44 
224 SLOPE=F(DIFF-DIFFO)/ (RWC J) -RWO) 
RWNEW=RWCJI-DIFF/SSLOPE 
RADOT(CJ)=(RWNEW-RWREF )/DELTAU 
VAP (J) =QUID(J+1)-QUID(U)-RWDOT( J) 
Dir FO=DIFF 
RAOSRW( Y) 
RwW(JSRWNEX 
IF (LTIMES+ 3g )226,44,44 
226 PRINT 227, 0UNBAL(J),DERIVHS J) 
227 FORMAT(LSHRESOILER LOOP,2F20.8) 
SO 7G 30 
22> DP) (226 [= 
KtivJ=ie 
2260 X€Cladev=TAC lid) 
255 NIALLY=NTALLY <1 
IF ONTACLY“NPRINT) 4,217,127 
17 PRINT US»NIALLY,DELTAU 
18 FORMAT(CLZ7ZAHLINCREMENT NO. = I5/34HOINDEPENDENT VARIABLE INCREMENT = 
a ae erie “4 
PR2RiNt 29 
19 FORMAT(L16HOSTG MOL FRA SUM H UNBALANCE H DERIVATIVE TEMPER 


i-NC 
Ce oO eee Db 
( 


LATURE LIQUID FLOW LIQ ENTHALPY VAPOR FLOW VAP ENTHALPY ) 
SS TS=6 
JLIM=ENS ed 


ae 27 Je2.sciM 
PRINT 25,45576,SUM(J),QUNBAL (J) ,DERIVH( J), TEMP (4),QUID( JU) ,QUIMOH(Y) 
1L,VAP(J),VAPMOH( J) 
2G FORMATCI4,6E145.8) 
o7 ot ice JSiGs] 
oe 272405,9R 
“On, CONJENSER DUTY = €29.6/27HQREBOILER DuTY = E45.8) 


5 


as 


274 


may 


wma ar Odo) CW 
oe 


-«t Ute 


tg} ge 


SiGe (XC1,J),151,N0),RW( 4) -RWDOT (J) 
5 STAGE NO. = 13/(10F15.10)) 


“7 


“1.0m 30 
~ Sh Cu 
~ 
. 
C 
a 


Cc 
wt 


29 


“PU ta 

a Ns 

G) < if NIG) i 
3- 


: 


Mee aIN 0 eNO ree meso 3 


LPCNTACLY=NONTROW) 45 29,50 


Sho Soe, tose 
END 
c 
G SUBROUIINE FOR BUBBLE POINTS 
ie 
SUSROVTINE 3L8PT(T, i, GPR 
DIiMenSi ON GoOpokG iGo Al io oc eos cbr) 
COMMON 9G. LUA ses oat 
COULLKFE (AK 32,0,7) SEXP CAS (T4555 .56549+C8(7T+469.09)) 
KTIMES=2 
5S SUMY=0 <9 
DOCS eee 
YeEQUL En Cat, 9504) C10) 5) peu lune) 
4 SUMY=SUNY+Y 
iF CABSF (SUMY-2.0)-BPERRIE,0,5 
5 KTIMES=KTIMES-2 
ip (KL LMeES) 7.656 
6 SUMYO=SuUMY 
TO=T 
i eer) 
GO Pou 
7 SLOPE=(SUHY-SUMYO)/(T-7~°) 
TN=(C2.6-SUMY)/SLOZME)+T 
SUMYO=SUMY 
Loe 
T=IthN 
Go, 107s 
8 REIURIN 
oND 
END 
OC AND OR INITIAL VALUES 2650, 08752567 135687.02693748 
Gr AT CONVERGENCE 43247.25945198 
ALPH SET Bigeye 1 Ce a ~SSO 226 
nNOLDUS RATIOS 
006000 20.00c000 ~9859G67 .994619 .999095 1.000669 .84492190 


OROGRAM DYNAMIC 
C THIS #290GRAM HAS BOTH CONDFNSER 24ND REBOILER DELAYS AND ACCUMULATOR 
OTMENSTAN X(705 90). TKC10,50..X1F¢ se S00. YIE (10.90). XTSSF(10.50).7! 
4SSF(10,50),ENTHK(10),FNC(10).DLYV(50).,0LYHV(50),DLYHL(50), 
? ENTHI (10),ENTHV°10),FNTHW(10),FK(10),ALPHA(S),BETA(S),GAMMA(S), 
3NFERIVM(10),TEMP(50).VAP(59),QUIN(50),FL(50).FV(50),SSFL(50),SSFV(5 
4n),SSPL(5n),.SSPV(50).SUM(50).Y¥A2MNH(50),QUIMOH(50),DERIVH(50),QUNB 
SA! (50),HFL(50),HFV(50),HSSFL (39), HSSFV(50),RW(50),RWDOT(50), 
6RHO9ND(1H).CON(1N). FFEN(10),QUIDK(12).4(10),B8(10),C(10)-DLYX(10,.50) 
7, NELAYL(50).DELAYI H(50).DELAYX(10.50),P(20,40) 
COMMON AUIDX,A,B,C 
EQUILKF(A,8,C.1T) =EXPFCA/(T+460.0)48+C#¥(T+460.0)) 
THALPFCY,Z,T)=¥eToZ 
RFAD 2,NC,NS,JFEEN,.NCNTRO! 
2 FORMAT(415) 
READ 333, (ALPHA(1),BETA(1),GAMMA(I1),121,3) 
333 FORMAT(AF12.6) 
1 FORMAT(AF10.4) 
JILIMSNS+2 
BEAD tet (X01. J) cleat nC). Jsi1,Je1M) 
READ 2sCUXTF C1. J), Let NC). Je2.NS) 
SEAN fe( (VIF C13). 1St.NC).J=20ns) 
READ te CXISSE (1s )s Fe NC lo d= 20NS) 
READ ZeC(YISSE( 1, J)sT22,NCV,Je20NS) 
READ 1,.¢TEMP(J).J=1-JL IM) 
READ 1,¢(VAP(J),J24,J5LIM) 
READ 1,¢QUID(J),J=1-JLIM) 
RFAD 1,°FL(J),J22.NS) 
RFAND 1, ¢FV(J),J=22.NS) 
RFAD 41,(SSFL(J),J=2-NS) 
READ 1,¢SSFV(J),J=2.NS) 
RFAD 1,(SSPL(J),J=2-NS) 
RFAN 1,(SSPV(JU),J=2-NS) 
READ 1, ¢(VAPMOH(J).J=1.JI IM) 
RFAN 1, (QUIMOH(J).J21.JI IM) 
RFAD 1,¢(HFL(J),5=9,NS) 
READ 1, (HFV(J),J322,NS) 
QEAD 1. (HSSFL(J).U22,NS) 
REFAN 1, ¢(HSSFV(J).J=2,NS) 
RFAN 1, (RWC), J=2. JIM) 
RFANL» (RWNOT(J),J=2.NS) 
READ 1, (ENTHK(I),121,NC) 
REAN 1,.¢ENTHL(I),121,NC) 
READ 1, ¢(ENTHV(!1),121,.NC) 
READ 1, (ENTHW(I),7T21.NC) 
READ 1,(0AC(1),121,NC) 
READ t,(hC lie lets NC) 
RFAND 1,(C(]),121.,NC) 
RFAN 1,NELTAU,BPERR.DFLFRR 
READ 1,(RHOQD(I).12=1.NC) 


QC = 13297.23948/93 


DO 560 121,50 (Wivee CearOrssec/<$ Fert 
DFLAYL (1 )SQOUID(3) DELAY INTECUAL (£0 5%‘) 


DELAYLH(L) =QUIMOH(3) fot VAMR FLOW TO COMOENER- 
9 


561 
564 


100 


445 


Aas 


44 


ay 


Wh Yh a oe tees 
Di YaVvVeLysVAPMIHCNS ) 
Mh VaLCLVS I POA CNS «2 ye 
04 3564 124,50 ~ 
Pol ae Tell Pek Cl se 
Dee Jax Cl eee ie 
CONT INIE 
ene 4.4? 
= eS 
AREAS 25..49 
FN=0.27426 
Ci =0.42 
GR=2s9.2 
RAD Sree) 
OY Teo Slaten 
RHOLERANL + XCT,JFEFD) *RHNAD(I) 
CONSTHCT¥SGRTFL2.«GR)«eAFIR 
NENIMETCONS: *f HSI 
MTSHW + (TD/NENOM) w%, 6667e(GUID(CUFEED) )**.6667 
WI RSSSRHOL wAREAuHT 
Al PHZARFAwRANSELRSS 
BETSAREASWLRSS~(FA/SCONST ) ww. 6667 
1 ASQEI SRHAL PR BFT 
PQINT 6.LAREI.,ALPH,RET 
FORMAT(A6,9F12.6) 
NN 9 J=ER,LNS 
RHOL =9.90 
io. Si a 
QHOLZERANI. + XCI1,J)*RHCON(T) 
RWC JS=Al PH¥RHOL + BET (RHOL*=.5333)*(QUID( IU) **, 6667) 
PRINT TACC, «RAC Je Jal. J) tM) 
PARMAT(15H HOLDUP RATIOS ,/,12F10.6) 
SHi = J0.1"peELTAU 
NTALLY=N 
NPRINT=A 
JPATH=L 
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220 ‘LVASY=5n | 


U 
69 500 1=1,49 DELAY (I 7€28 VBL a, 





DI YVCLVARY)EDLYVCI VARY-1) 
ni YHV(LVARY)=DLYHV CI VARY-1) lal COMPEN SER VAPOR CME. 
Di YHLCLVARY)=DLYHI CILVARY=1) ; 
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SUBROUTINE FOR RURBIE POINTS 


SUBROUTINE RUBPT(T,L,BPFERR) 

DIMENSION QUIDX(29),A(10),8¢19},06(10) 
COMMON QUIDX,A,B,C 

FOUTLKF (A,B,C, T)=FXPFCA/(T+460.0)+3+C#(T+469.9)) 
KTIMFS=14 

SUMY=0.n 

NO 4° T=4 51 
YeCOUILKFCAC(),8¢1) ,6C(15,7 )eQ0l0Kt 1) 
SUMY=SUMY+*Y 

IF CABSF(SUMY-1.0)-BPERR)8.8,5 
KTIMFS=KTIMFESe1 

TF CKTIMES) 7,6,6 

SUMYOSSUMY 

Ths tk 

ToT+19.N 

GO TO 3 

St OP FHZCSUMY-SUMYO)/(T-TO) 
TN=(0(01.N-SUMY)/SSLOPF)+T 

S'IMYO=SUMY 

TN2T 

T3TN 

Ga 10 Ss 

RETURN 

END 

END 

AR INITTAL VALUES 11849.03732467 13687.02693748 
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APPENOIX IE 
COLUMN PARAMETERS FOR HYORAULIC EQNS. : 





FEED RATE = 250 gal/min a 
COLUMN DIAM. = 6 FT. 
ACTIVE TRAY AREA = 76% TOWER AREA 


WEIR LENGTH = 77% OF TOWER NG 


WEIR HEIGHT = 0.25 FT 


CALCULATED VALUES 


‘= 4.62 FT WEIR LENGTH - ig 
A= 21.49 FT? ACTIVE AREA | =e - 
Y.# 0-25 WEIR HEIGHT | | 

Fs 16.576 MOLE/MIN 

C= 0.42 WEIR EQN CONSTANT 


THESE VALUES ARE THEN USED TO CALCULATE COLUMN CONSTANTS 


| a; 2; | 
: 3 3 LIQUID HEIGHT 
a (7a Tei (,) ON STAGE; 


(Wirdss = GAYS HOLDUP REFERENCE STAGE 
| (FEED STAGE USED) 


cz. ni ae 
 UWerlss 


i 
nt. (Winss © g77 


}constanrs FOR EQN |2 
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APPENDIX IV 


COMPUTER FLOW DIAGRAM FOR NUMERICAL SOLUTION OF GENERAL 
STAGE EQUATIONS 


Summary of equations is as follows: 


Kaj ¥45/% 44° explA, / T, +B, +C,T 3) 3) 
ot eae or alee (*) 
iy Ry, {(Sinput - output) - Ra Xi (7) 
@ 1 e 

h, = a ((Sinput - Hout put) aveey - Ryjh4/ (8) 
Ry = ies + BF, LS (9) 
Raj = (Linput - Youtput) ass aeraas (10) 
h, = a,T,. +b (12) 


a a | 


By Pky (15) 


j r(a,T, ++ ae + are KK ac) (20) 
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APPENDIX V 


FORTRAN 60 PROGRAM OF LEAST SQUARES FIT OF MODEL RESPONSE 
AND A SYNTHETIC TRANSFER FUNCTION 


First program is for a two parameter fit of transfer function 


-.338 
e 


where G(s) is the transfer function of C, Component of overhead 


4 


product for a step in C, feed composition. 


4 
Second program is for a three parameter fit of same transfer 


function with an arbitrary gain adjustment, K, i.e. 


Ke 8 
Gs).= (7,841) (7,541) 
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PRUGRAM Lsis -CrWwo PReaAMEeTE Fir ) 


WRITE FUNCTION SUBPROGRAM 79 CALCULATE THE FUNCTION 10 8E LEAST SQUARED. 
A IS DIMENSIONED 10 AND CONTAINS THE PAHAMETERS 10 BE VARIED SUCH AS TO 
MINIMIZE THE SUM VF ERRORS SQUARE x, XB, AND XC ARE THE IHREE INDEPENDENT 
PAI TERS A YOU MAY WISH 10 USE ONUY x) NOR IS THE FUNCT ION ROMBCoEm USED FOR 
BRANCHING TO DIFFERENT PARTS OF THE EQN SJBPRUGRAM WHEN SEVERAL FUNCTIONS 
ARE TO BE FIi (FG- SEVERAL JOBS TD BE DONE). 
THE DIMENSIONED ARRAYS HAVE THE FOLLOWING MEANING-- 
A(1Q) CONIAINS THE PARAMETERS wHICH ARE TO BE VARIED 
X(200)+ XB(200), AND XC(200) CONTAIN THE INDEPENDENT PARAMETERS. 
Z(200) CONTAINS THE OBSERVED VALUES OF ‘HE FUNCTION (THE OBSERVED DEPENDENT 
PARAMETERS). 
FINCR(i0) CONTAINS THE MAGNITUDE SF THE INCREMENTS FOR THE PARAMETERS A(10) 
SO THAT (HE PRUGRAM MAY TAKE NUMESICAL DERIVATIVES WITH HEASONABLE ACCURACY. 
R(200) CONTAINS THE DIFFEREVCE BEIWEEN OBSERVED (Z7(200)) AND CALCULATED 
E(10) CONTAINS THE ESTIMATES OF THE ERAOYS OF THE PARAME/ERS A(19) AFTER 
ITERATION, 
INPUT 
FIRST CARD - WILL BE REPRODJCED Ov OUTHUI USED FOR LABELING. 
SECOND CARD FUKMA (12,13,12,13,E10.2) 
[R= NUMBER OF PARAMETERS TO BE VA-IED 
1S NUMBER OF POINTS 
NOF = FUNCTION NUMBER (3EE ABOVE) 
NINP = NUMBER OF INDEPENDENT PARAMETERS 
EPSIL = IS USED AS A CRITERION FOR CONVERGENCE. IF THE RELATIVE VALUE 
OF THE RESIDUAL CHANGES BY LESS THAN EPSIN IN TWO SUCCESSIVE ITERATIONS, 
CONVERGENCE IS REACHED. 
THIRD CARD(S) - FURMAT (5E14.6) - CONTAIN YOUR ESTIMATE OF THE PARAMETERS, 
A(10) 
FOURTH CARD(S) - FORMAT (9614.6) - CONTAIN THE INCREMENTS FOR THE A(10) FOR 
OBTAINING NUMEKICAL DERIVATIVE. 
FIFTH CARD(S) -FORMAT(5E14.6) - CONTAIN Z(1), OR THE OBSERVED. 
DEPENDENT VARIABLES. 
SIXTH CARD(S) = FCRMAT(SEL4.6) - CONTAIN X(1), OR THE VALUES OF THE 
FIRST INDEPENDENT VARIASLE. 
SEVENTH CARD(S) A::D EIGHT CARD(S) ARE READ IN ONLY IF NINP [TS GREATER THAN 
UNITY, AND CONTAIN XB(I) AND XC(I). THE SECOND AND THIRD INDEPENDENT VAR] ABLES 
SACKING OF JOBS IS PERMITTED - THE PROGRAM STOPS WHEN IT READS A ZERO FOR ip 
(SECOND CARD INPU'), THUS SEVERAL BLANK CARDS AFIER LAST yog Is 
SUFFICIENT TO STO+ ROUTINE. 
DIMENSION A(10), x(200), XB(200). xC(200), 2(200)» FINCR(10), 
18(200), £(10),(40),2R(10),W(10) 
1 READ 106 
READ 100, IK,IN,NOF.NINP,EPSIL 
IFCIR) 200-200-190 : 
10 READ 104. (ACL » | = 4,1) -. 1 foe 
READ 101+ (FINCR(]). I = 1,1R) SxPesimenrar pots XC) ie. xE-dogre) 


J ) i es 
READ 102. (wkiI}s]zisi8) mem Rees 0 SEVER Lereawons Detywomi@ee 
READ 102-(wWt1l).}21,18S) ae CRISTO MFRE Wweréwrs 
PRINT eotUCl erste tS, (RCL) sda, (50, CWUl)s Pais lS) 
2 FORMAT (16HU +R Ww INPUTS,./,(7E12.6>/5) ) 
READ 103- (X(T , 1 = 1513) deme YRewrs om SBP“ corr AM 
DO 18 [s1,I195 
WW = X(T) - 4.9 
ZEL ys 000 
DO 19 J= 1» IS 
19 21) ga Ps Ce aa) ee WW ea OK FB) NER THE x (s) vacugs 
18 CONTINUE 
IF CNINP-2) 20.,15,15 
Lo -WEAD 1042. AXE CIS5ac = 2,985 


FO 


16 
20 


25 


26 
30 


100 
101 
102 
103 
104 
105 
106 


107 


108 
109 
110 
111 
112 


IF CNIN 
READ 1 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 


P-3) 20.16.16 

05. (xCcid, | = 2,18) 
106 

109 

107 

108, (IR.IS,NOF,.NINP.EPSIL) 
109 

110 

109 

Lil, (ACI). | = 2,9R) 
109 

112 

109 

Ji1l, (FIWCRCI). 1 = 1,IR 
109 

1135 

109 

114, (Z(id. IT = 1.18) 
109 

115 

109 

116.5 (Xi), I 
109 


1,1S) 


IFCNINP?2) 30,25,25 


PRINT 
PRINT 
IF CNIN 
PRINT 
PRINT 
CALL L 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 


117, (xB 1), | = 1.99) 
109 

P-3) 30,¢6.26 

118. (xCil), J = 1-95) 
109 

es ne cee ERE NODES Aone Ine 
119 

109 

1205. (ACi)>. | = 1.TR) 
109 

124 

109 

120, (€cl), T 2 1,7R) 
109 

121, NOI'R. RRQ? 

109 

122 

109 

125- (R(id, | = 1,198) 


Go TO 1 


FORMAT 
FORMAT 
FORMAT 
FORMA 
FORMAT 
FORMAT 
FORMA: 


FORMAT 


FORMAT 
FORMA! 
FORMA: 
FORMA! 
FORMAT 


1) 


115 FORMAT 


(12.-13,!2,13,.610.2) 
(5614.6) 

(5E14.6° 

(5614.6, 

(5614.6) 

(5614.6° 

(80h 


(99h NUMBER Or PARAMETERS NUMBER UF POINTS FUNCTION 


1 NUMBER NUMEFERQ OF INDEP PARAM E®STLON) 


dee coe ea ed Pees Re rest Mee 

H ) 

(26h THe PARAMETERS FED IN ARE) 

(1H »5614,6) 

(55H THE INCREMENTS FOR OBTAINING NUMERICAL DENIVATIVES ARE 


(34h THe OBSERVED VALUES TO BE FIT ARE) 


¥/ 


QQ 


AagANAIgaAAaAN 


114 FORMS? (5614.6) 

115 FORMAT (31H THE INDEPENDENT PARAMETERS ARE) 

116 FORMAT (1H .56€14,6) 

117 FORMAT (1H ,5€14,6) 

118 FORMA] (1H ».56€14,6) 

119 FORMAT (38H THE BEST VALYES OF THE PARAMETERS ARE) 

120 FORMAT (1H ».7£17.10) 

121 FORMAT (24H NUMBER OF ITERATIONS = ,12,50H 'HE SUM OF THE SQUAR 
1—S OF THE ERROFS IS NOW = ,£12.6) 

122 FORMAT (53H OBSERVED FMINYS CALCULATED VALUES OF THE BEST FIT ARE) 

123 FORMAT (1H .9612.5) 

124 FORMAT (66h ES'IMATES OF THE ERROR IN EACH PARAMETER ARE (STANDARD 
1 DEVIATION)? 

200 STOP 
END 


SUBROUTINE LEAST(IR,IS.A,X,XB,XC,Z,FINCR,EPSIL»NOITR»RRO,NOF,RDE) 
IR = NO. OF PARAMETERS, IS = NO, OF POINTS, A IS ARRAY OF PARAMETERS, 
X IS INDEPENDENT VARTABLE( Z DEPENDENT. FINCR IS ARRAY OF INCREMENTS 
PARAMETERS. EFSIL IS -FRACTIONAL- ERROR CRITERION. NOITR IS NO. OF 
ITERATIONS REQUIRED (UP [C 10). RYO = SUM OF SQUARES OF RESIDUALS) 
NOF IS NUMBER OF HE FUNCTION USED IN -EUN-. 
E IS THE ARRAY OF ESTIMATED ERRORS IN THE PARAMETERS 
DIMENSION AC102, X(200), XB(200), XC(200), 2(200). FINCR(10), R¢20 
10)» D(200,10), DT(10-200), DEL(10), DS¢10), DPI(10,10), DP(10.10), 
2€(10) 
NOITR = 
0020 (o=ad7 hs 
20 RCI) = ZCI) - EQNCA,X(1),XBC 1). XCC1),NOF) 
PRINT 4, (K( 1). [=1,18) 
25 IF (NOITR-9) 130,130,101 
130 NOITR = NOIIR + 1 
DO 220 J = 1, TR 
ACJ) = ACJ) *FINCROY) 
OG 15) ob = 440s 
15 OC}eJ) 2 EGNCA.XCT),XBC1).XC(L)»NOF) 
ACJ) = ACU d-FINCR(Y) 
220 CONTINUE 


0 


DO 30 J = 1-158 
CONST = EQN(A,<(7).xXB(1),XC(1).NOF) 
DO 30 J = 4-IR 


30 DCI+J) =s (D(1,.5)-CONSTI/FINCR(Y) 
PRINT Se ClLDCl se) aJeio IT RY bet, 1S) 
3 FORMAT(29HDERIVATIVES OF EACH PARAMETER,/,(4F10.5,/,)) 


DO $5 [ = 4-1S 
DO 35 J = 1-IR 
SoeD teal) = PUP) 
DO 36 I = 1,5IR 
DO 36 J = 1,IR 
DP(I,J) = 0-0 
DO 36 K = 41;1§ 
36 DP( I,J) = DP (CT. J)+DTC I,K) eD(K,Y) 


CALL GAUSSS (1R,1.00E-30,DP,D-1,KER) 
IF (KER-1) 120,37,120 


37 DO 40 1 = 1,515 
DO 38 J = 151R 
DS(J) = 0.0 
DO 38 K = 1-IRF 
38 OS(J) = DSCI) +E PT (CU, K ODT (K. {) 
DO 39 L = 1-IR 
39 DI(Leld = DSCL> 


40 CONTINUE 


gz 


OF 


110 
19 
320 


120 
1001 
100 


140 
150 


11 
12 


13 


14 


15 
20 
30 
31 


32 


33 


DO 110 [| = 1. iR 

DEL( I) = 0.0 

00 110 J =i. '§ 

DELCI) = DEL( LT +DTCI»c ety) 


DO 10 J = 1-1 

ACI) = aCl)*DELCT) 

00320 I = 1-18 

RCI) = Z(t) - EQN(A. KOT), XBCE)-XCCL). NOP) 
PRINT 4,¢R(1),121,15) 

FORMAT (16HDOIFF IN OS8SS*CALC./,7F10.5) 
RRQ = 0.0 

DO 50 I = 1-18 

RRO = RRG+KI 1) e#2 

CRES = ABSFORRO=-RRP) - EASTL e RRO 
RRP = RRQ 

IF (CRES) 100,100.25 

FORMA: (20h CONVERGENCE FAILURE) 
PRINT 101 

GO TO 100 

PRINT 1001 

FORMAL (16H MA’RIX SINGULAR) 
FPS? = [3-1 

DO 150 [| = 1,1T* 

SUM = 0.0 

DO 140 J = 1-J> 

SUM = SUM*D(JU.:)ee22 

ECT) = SORTFCR Q/(SU4eF I S51)) 
RETURN 

END 


SUBROUTINE GAUSSS(N,EP. A,X. KER) 
DIMENSTON A:10.10). (10,10) 
DO 1 [#1.N 

DO 1 J21-N 

X{I2J)29.0 

DO 2 K=1.N 

xX(K2K) 214.0 

DO 34 Lei.h 

KP=0 

220.0 

DO 12 Kel-h 

IF CZ-ABSF CACK.1 ))111-12,12 
Z2ABSF(CACK.L)) 

KP=K 

CONTINUE 

JF (L-KP)13,-20.20 

DO 14 J=L>N 

Z2eAClL.J) 

A(L»~J)SACKF.J) 

ACK®,J)sZ 

DO 15 J31.N 

Z2eX(LoJ) 

K(L>/J)SXCKF.J) 

M(KP, J) 22 

IF CABSF(CAC(L»L))-EP)50.50.30 
JF (L “ND 31534,34 

LP1=L+1 

DO 36 K=LP1.>N 

[FC ACK,LI)$2,36,32 
RATIO=A(K.LIZAILoL) 

DO 33 J=Lh1.N 

ACK, JI ZACK, IJ) -FATIQ*A(L, J) 


00030 


00030 
00040 
00050 
00060 
00070 
00080 
00090 
00100 
00110 
00120 
00130 
00140 
00150 
00160 
00170 
00180 
00190 
00200 
00210 
00220 
00230 
00240 
00250 
00260 
00270 
00280 
00290 
00300 
00310 
00320 


DO 35 Jzei,N 00330 


35 X(KeJ)SX( Ke) -KATIO#R®X(L, J) 00340 
36 CONTINUE 00350 
34 CONTINUE 00360 
40 DO 43 [=1,"% 00370 
LI=Ne1-] 00380 
DO 43 J=i.N 00390 
IF CTI-N) 41,493,453 00410 
41 JIPL=AlIo1 00420 
DO $2 KsliFi,N 00430 
42 SsSe*AClTI,Ki*x(K,J) 00440 
OS MCT IIe Cx Clad yeSI7k( 15.14) 00450 
KER=1 00460 
REI!URN 00470 
50 KER=2 00480 
RETURN 
END 00490 


FUNCTION EGN(A.X,XB-XC>NIF) 
DIMENSION A(10 





EQN=.052*EXPF(-S.eXd/C (ACL) OX%1.) e(A(2)eXOl. exey 
FUNCTION TOBE CIT 






U RR WwW INPUTS 
-S21390E-01 .189860E-01 .940000E-02 .463800E-02 .148000E-02 .530000E-04 -O00000E+ 


°254460E-01 .1292346+00 .297077E+00 .500000E+00 .7029226+00 .870766E+00 -974554E+ 


-647420E-01 .1359853E+00 .1909156+00 .208980E+00 .1909156+00 .139853E+00 -647420E- 


ef 


SYN UF XKFR FCI UF Fi wm Keer (o35)/05°Se1)(Gn01) 10 . UP PROD KES® 16 FEED PERT 


NUMBER OF PARAMETERS NJMBEL OF POINTS FUNCTION NUMBER NUMBER O» INDEP 
2 7 1 1 


1HE PAKAMETERS FEU J ARE 
,$67000E*03 41000002 

THE INCREMENTS FUR OBTAINING NJMERICAL DEnIvailVES AnE 
.100000E*01 -10¢000€+01 

THE OBSERVED VALLES O Be “I! ARE 


-7789B81E-02 -161892E-02 -591618E-03 -277002E-03 -148540E-03 
-8673645-04 -537U52E-04 


ThE INDEPENDEN! FARAPETERS ARE 


-100000E+01 -200000E+01 ~-300000E+01 -400000€-+01 -5S00000E*01 
-600000E°01 -700000E+01 


DIFF IN OBS-CALC 


-30018 ~000u1 .00000 -.00000 -.00000 -.00000 -.00000 
DERIVATIVES OF EACH PARAMETES 
-.00002 -.00C05 -.00000 -.00002 
-.00000 -, 00001 -~.00000 -.00000 
-.90000 -.00000 -.00U0U -.00000 
-.00000 -.00060 
DIFF IN OBS-CALC 
-.00001 -,00001 .00000 .00000 -.00000 -.00000 -.00000 
DERJTVATIVES OF EACH PARAMEIER 
-.UV0002 -.000uU5 -.00000 -.00002 
-.00000 -,00001 -.00000 -.00000 
-.00000 -.00000 -.00000 -.00000 
-.00000 -.00000 
DIFF IN OBS-CALC 
-10000 -.00000 .00000 -00000 -.00000 -.00000 -.00000 
DERIVATIVES OF EACH PARAMETER 
-.00002 -.00005 -.00000 -.00002 
-.00000 -.00001 -.00000 -.00000 
-.90000 -.00600 -,.00V000 -.00000 
-.00000 -.00C00 
DIFF IN OBS-CALC 
- 30000 -.60000 .00000 -00000 -.00000 -.00000 -.00000 





‘RE SEST VALUES Ch THE PARAMETERS ARE Two patarmeres piy 
$441729345E+03 .45847984888-02 ou ey 
SS a 


ESTIMATES OF THE ERAIR JN EACH PARAMETER ARE ‘STANDARD DEVIATION) 
of VieNMTI HEPES © Od9e rs Ete! 


- 


59 


NUMBER OF TiERATIUNS = 8 THE SUM OF THE SUJARES OF THE cRRORS IS NOW = .289€ 


OBSERVED MINUS CALCULATED VALUES OF !HE BES: FIT ARE 


-86381E°07 -.15422'-05 .27661F-05 .14834E-05 -.77337E-06 -.24053E-05 -. 330475 
STuP 


TIME, 2 MINUTES AND 122 SeCONDS 


56 


PRUGRAM LS S =m (3 PARAMETER eir) 


WrITE FUNCTICK SUbPROGRAM 1) CALC ILATE BE FUNC JUN 10 BSc LeAS SQUARED. 
A TS DIMENSTUNED 10 AND CUNTAINS HE PAW-AMETEXS 1G BE VanlED SUCH aS "9 
MINIMIZE THE SUM  F ERRJnS SQUARED ©, 48. AND KC ARE THE [ower INDEPENDENT 
PARAMETERS (Yl. May wISd TO USE ONLY ©) VO TS .He FUNCTIUN NUMBER © USED Fi& 
BYANCHING TO DIFFEWENT 2ARTS OF THE EQN SUdDPRUGTAM WHEN seveda! FUNCTIONS 
APE TO BE FI! (FO. SEVERAL JOBS 19 BE DANED. 
THE DIMENSTONEU Ar RAYS HAVE THE FOLLOWING “MEANING- = 

A(10) CUNIAINS THE PARAMETERS -HICH ARE TO BE VARIED 
X(200). XB¢200U),. 4ND xXC(200) CONTAIN Tet INDEX ENDENT PAWAMEI EWS, 
Z2(200) CONTAINS oe OBSESVE) VALUES OF 'HE FUNCTION (THE UBSenwVED DEPHNDENT 
PARAMETERS). 
FINCR(10) CONTAIN: THE WAGNITUDE OF THE JNCREMENIS FGR THe PARAMETERS A(10) 
SO THAT (HE PFUGHAM MAY TAKS NUMEXTICAL DERIVATIVES WITH wEASONABLE ACCURACY, 
R(200) CONTAINS ttt DIFFEREVCE SE: WEEN UBiscRVED (2¢200)) AND CALCULATED 
E(10) CONTAINS THRE ESTIMATES OF THE ERRORS OF THE PARAMEIENS A(10) AFIER 
TTERATION. 

INPUT 
FIRST CARD - wiLL BE REPRUDICED ON OUTRi: VSED -OR LABELING. 
SECOND CARD FURMA'([2,)5.12.133,E10.2) 
IRs NUMBER OF PARAMETERS [TO BE VAXIED 
Is= NUMBER OF *OI::TS 
NOt s FUNCTION NUPFBER (35EE ABUVE) 
NINP = NUMBER CF [INDEPENDENT PARAMETERS 
ePSIL = IS UStD A~ A CRITERION FO- CONVERGENCE. JF THE RELATIVE VALUE 
Ur THE RESIDUAL CrANGES BY LESS THAN EPSIN IN Twu SUCCESSIVE JTI!ERATIONS. 
CUNVERGENCE IS REACHED. 
1HIQD CARD( >) - Fi RMAT (5614.6) = CONTAIN *OURK ESTIMATE Ub THE PARAMETERS, 
A:.10) 
FIRTH CARD(S) - FORMA! (9E14.6) - CONIAIN THe INCREMENTS FOR IHE AC10) FOR 
OBIAINING NUMEXITCAL DERIVATIVE. 
FIFTH CARKD(S) -FO'MAT(5514.9) = CUNTAIN 201), O% IME OBSe VED, 
DEPENDEN! VARIABLES. 
SIXTH CARD(S) - Fi wMAT( 9614.6) - CONTAIN Kid). OR THE VALUES OF THE 
FIRST INDEPENDEN: vARIAQLE. 
SEVENTH CaRD(S) AV‘D EFIGHE CARD(S) ARE xEAD IN ONLY ITF NJ WR IS GREATER THAN 


UNITY, AND CONIAIN x¥BCT) AND XCC J). THE SECUND AND TWIRD INDEPENDENT VART ABLES 


STACKING OF JOBS 1S PERMITTED - THE PROGRAM STOPS WHEN |]! READS A ZERO FOR JR 
(SECOND CARI INP 9, THIS SEVERAL BLANK CaHDS AFTER L&aS!T JOB IS 
SuPFICIENT 'O STuU- ROUTINE. 
DIMENSION &4(10>, x¢€200). x8(200). x¥C'200), 2(200), FINCR(10), 
1R(200), Eci1U). (10),%81°0).w(10) 
1 READ 106 
READ 100+ InH,1~.NOF.VINP.ERSIL 
IF CIR) 200,260-10 
19 READ 101- (ACT. — = 1-12) 
READ 102, (hiwere cy). | = 1,1) 
READ 102. Cu] elas 1S) 
READ 102. (nf&(:).fP21-I3) 
READ 102-¢wt]). 1221-135) 
PRINT 2,¢UC1),.,214,18).C(RVUCT)-T=1,7S)-Catt).y=1.-198) 
2 FORMAT (16h -R_ A INP ITS e/4-(7E12-.6+/2) ) 
SEAD 103. (ACT . | = lel) 
N0 18 [=1.1> 
wwo2 X€1) = 4.6 


Z(1) = 9.0 
dO 19 Jz 1. IS 
1¥ ZO) = Jibs © © CJ CUCL IORI J) eean 


18 CONTINUE 
IFCNINP=2) 20,15,19 
15 wEAD 104. (‘aRt )d, | = 1-19) 


4 


16 
20 


25 


26 


30 


100 
101 
Pi2 
103 
104 
103 
106 
a 
107 
il 
108 
109 
116 
111 
rie 
bl 
13 


IFONINP]-3) 20.16.16 
REAO 105. A Gels I = taks:) 
PRINT 106 
PRINT 4109 
PRINT 107 
PRINT 100+ (IK-I1S,NGr ,NIVP,EPSIL) 
PRINT 109 
PRINT 110 
PRINT 109 
PRINT Tite WAC) se ol TR) 
PRINT 109 
PRINT 112 
PRINT 109 
BOUNTY ot tdci Ceol Moy seers al Re 
PRINT 109 
PRINE 113 
PRINT 409 
PRINT 224, “20i)y [| = Ls7S) 
PRINT 109 
PRINT 115 
PRINT 14109 
PRINT Tios Ci) b= oa lSo 
PRINT 109 
IFONINPe2?) $0,.25,25 
PRINT Divs (X84 195 1 = dels) 
PRINT 109 
IF(NINPS3) 30,26,26 
PRINT $2265 txCt])}, 4] = 1415) 
PRINT 409 
CALL LEASTCIF,'S,A.X.xXB,XC,Z,F INCR,EPSIi,NOLIR,RRG.NOF RE) 
PRINT 109 
PRINT 119 
PRINT 109 
PRINT 2120- tACl). T = 1,78) 
PRINT 109 
PRINT 124 
PRINT 109 
PRINT 120, ‘ECi),» I = 1,7R) 
PRINT 109 
PRINT 2215 wWOI R,. RI 
PRINT 109 
PRINT 122 
PRINT 109 
PRINT 123. (Ri), [| e5ie FS) 
eG Fa aie | 
FORMA] (12,13.12,13,¢10.2) 
FORMA! (5614.6. 
FORMAT (5614.6) 
FORMA! (5614.6. 
FORMAI (5614.6 
FORMAT (5614.6 
FORMAT (80K 
) 


FORMAI (99F WUMBER OF PARAMETERS NJIMBER OF POINIS FUNCTION 


NUMBER NUMBER OF JNDRP PARAM E25T7LON) 
FORMA: (84-12,¢4K,13.17%-12,19%512,12%,+10.2) 
FORMAL (1H » 

FORMAT (26h THr PARAMETERS FED IN AXE) 
FORMA! (1H + 5t14.6) 


FORMA: (55r THE INCREMENTS FOR OBTAINING NUMERICAL GERIVATIVES ARE 


) 
FORMA, (34h (He OBSERVED VALUES TO Se FIT ARE) 


$ 8 


a 


114 FONMAS (5614.6 

115 FORMA: (Sib “he INDE PE VCZN Phe AMET SH Ave. 

116 FORMA: (ih .5ri 4.6) 

117 FORMA: (im .56-24.6) 

115 FORMA! (1h »-56:4.6) 

119 FORMAT (36 TR Best vAL JES JF Tak -Ands-Tews ARE) 

120 FOKMA: (16 .7F17,.10) 

121 FOWMA: (24h AurwER UF i :etalloiua = .i2.,704 Me SuM oh tee S uae 
LES OF THe EFS5-S TS Vin = 1612.6) 

122 FORMA: (53m O68 ERVED VINIS CaLCuLatei vai ueS tb THE ge oil Fi Aw t) 

1235 FQORKMAI (1H .9F!2.5) 

124 FORMA! (66r ES JMATES Je THE E-RIQ J FOTH CA- AMET EX awE (> aN DARD 
1 DEVIATIUND» 

20) STQ° 
END 


SUBWOUTI ve Lea TCE Se AK eI AC eZee FC EMPSIL NUT TH NU NDt mE) 
I}? = NOe OF HPAW~ AME TERS, To = WO, 1 PQ Se A JT ARRAY Qh RPANWAMETERS, 
X TS INDEPE NDE’) ' ARJABLEC 7 DECE-Den:. -ivCwr 1S ARRAY Ob TNC: FMENTS OF 
PAHAMETERS, EFSTE TS FR ACTIONAL- E899 CAPE UN,. NOliw TS NO. OF 
PVeRATIONS ReGt bel (UR TE Ld), RQ = Sty Ul SQ ARES OF WESTDUALS) 
NOP TS NoMBew vb HWE FUNCTION USED IN -eEQN-. 
E 1S THE ARRA1 wb ESTIMATED =RRORS IN Ter PARAME TENS 
OIMeENSIO?# 4410 . x(200),. KB( 200). ¥5°200), 21200), Fi .CR(10), (20 
1N)d)- Di2ou-16). DVC1H-200)., DEL 10)- A810). De l(10,10). DPlCiO.190), 
2&(10) 
NOITR = 0 
no 20 | 1°TS 
20 RCE) = ZiT) - FUNGAL, COL). XBOCT)-XCCLI NOD 
PRINT @, Ce vd). Fet,tS) 
25 IF (NDITR-G9» 140.130.1201 
130 NOJTR = NU] ‘kK: 
DO 220 v= 1. th 
ACJ) = AtulrbhTE CQO) 
m DO 15 1 = 2-15 
15 Diode & EUNCA- CT), KBCL 9.00019, NOFD 
ACS) = ACs CROS) 
220 CONTINUE 
DO 30 1 - leis 
CONST =F curv A, CLV xdCLI-XCCLI.NOF) 
DO 30 J = fal* 
30 O¢CL»~J) = Cis d.. VOCINSDT ISK INCRCY) 
OQiInT. 3.6L vis. de Je 1 Reet. Ts) 
Z FORMA CA29MLER] .ATIVES JF EACH 2APAME'E RL 7, (4F-10.5,7/.)) 
BO'35 1 = 11S 
09 35 J: leh 
$a DOMCIe-]) = ube. ) 
DO 36 I - d-tr 
DO 36 J = 1-I* 
DP(I.U) - G.9 
DO 36 K = 11d 
36 DPCIT.~J) = Lert JGd*DTCLAdeD( KI) 
CALL Gaur. 3 ¢.6,1.008°3),0P.0 IKE) 
IF CKER-1) 120.357,120 


37 00 40 I J-I5 
bOI 38 J - Jein 
DSCJ) = O.6 
DO 38 K - 1-ik 
84 GSCI) = OSturep Pbtu. CdeDICK, 1) 
D0 39 L = l-in 
S99 DRCe el oes a OE 


4 CONTINUE 


10 


320 


120 
1001 
100 


140 
150 


ii 
i2 


15 


14 


13 
2) 
3o 
31 


es 


ee) 


DO 110 | = 1, 

DEL(I) = 0.0 

DDO? 110 J = 1 us 

DEL(1) 2 VEC +DIClecd ett J) 


DO 10 1 = 1aAlk 

ACI) = avcl)+Nei cl) 

RCId = ZC14 -=) EGNCAy (C1). XBCT) XC C1 32ND) 
PRINT 46, (R(1),7T21,18) 

FORMAT (16HDIFF IN O08S-C4_C,/,7F10.5) 
RRQ = Q.0 

DO S071 -=)1als 

RRQ = RRUski])+#2 

CRES = ABSFE(RRYU-RR®P) - EPSTL eR? 
RRP = RR 

IF(CRES) 100,100,25 

FORMAT (20h CU’-VERGEVCE FAILURE) 
PRINT 101 

GO TO 100 

PRINT 1001 

FORMAT (16H MA: RIX SINGULAR) 
FIsi-= ]$s-2 

DO 150 1] = 1-1 

SUM = 0.0 

DO 140 J = 1,15 

SUM = SUM+D\ 4,1 )e#2 

ECI1) = SukTFECRRFO/(SUMeF 1 31)) 
RETURN 

END 


SUBROUTINE GAUSSS(N.EP A,X. KER) 
DIMENSION 4:10.10), X(10,10) 


DO 1 121>N 

DO 1 J=1i-N 

X(1T,J)=0.0 

DO 2 K=1-+N 

X(KeK)=4.0 

DO 34 L=l1-hk 
KP=0 

220.0 

DO 12 KeloN 


IFC Z2-ABSFCACK,! 9))11-12,12 
Z=ABSFCA(K-L)) 

KP=K 

CONTINUE 

TPCLeKP I 1S. 20.20 

DO 14 Jer, 

L2ACE sD 

A(LesJIZACKFE> J) 

A(K®,J)22 

pO 15 Js]. 

Zax ed) 

X(LoJdISX( KFS) 

X(K9,J)2/ 

IFCABoF (ACL -L) -EP)50,90,30 
IF CL-NI31,949,3% 

LP1=L+1 

DO 36 Kelri-N 

IF CACK,L.)$2¢,56,32 

ad OSA UR AAD bl) 

vo - 83 Veber = 

ACK, JIZACKs VI -K ATI DAIL, J) 


Ge 


00010 


00030 
00040 
00050 
00060 
00070 
00080 
00090 
00100 
00110 
00120 
00130 
00140 
00150 
00160 
00170 
00180 
00190 
00200 
00210 
00220 
00230 
00240 
00250 
00260 
00270 
00280 
00290 
00300 
00310 
00320 


U 


35 


34 
40 


oF 
42 
43 


50 


RR 


DO 35 Jel. 
MCKsJIEX (ACH ATIOOK(L, J) 
CONTINUE 

CONTINUE 

00 43 [=1+4 

T1=Ne1-] 

00 43 J=1-h 

$30.0 

IF CJT I1-N)41-43,43 

PyPL1slle. 

DO 42 KzilFi.N 

GsseAC] ].Kiex(n,J) 

KOE SI SCecrt, »-SI/AcCITJ.,.] 1) 
KEF=1 

RETURN 

KEPs?2 

RETURN 

END 


FUNCTION ELwl Aen, XB>XCON IP) 
DIMENSION 4(10) 






END Custer ag Pir 

END 

® INPUTS 
.$22390E-01 .189660F-01 .940000E-02 46380UF-02 .148000E-02 
-254460E-01 612923546 :°00 .297077E+00 .500000e +10 .702922E+00 


~647420E-01 .139853E-00 .190915E-00 


208980E +0 


~190915E-00 


-530000E-04 
-870746E+00 


©139853E+00 


Leet 
CC S46 L 
CosSc 
00360C 
VOSZL 
OO<sel 
N0S3%C 
AQ40C 
v041U 
004206 
00430 
N044¢ 
06456 
90466 
90476 
u0480 


00490 


-OO0000E- 
~974554E> 


~-647420E- 


SYN OF XFR FC. UF Fi RM Keeer( -35)/01541)(85+1) 10 ‘OP FROD RESP 90 PEED PERG 


NUMBER OF PARAME'EFS NJIMBER OF PINTS FUNCTIUN NUMBER NUMBER OF INDEP 
3 fs 1 ik 


}HE PARAMETERS FtD [: ARE 
-344173E+03 -4956480ce+02 -520000E-01 

THE INCREMENTS FCR OCUTAINIVG NUMERICAL DERI VA'IVES ARE 
-100000E+01 -100000&+01 -100000E-02 

THE OBSERVED VALUES 9 BE rI! ARE 


-7789B81E-02 -161692E-02 -591618E-03 -277002E-03 -1498540E-03 
-867364E-04 -537052E-04 


THE INDEPENDEN! PARAMETERS ARE 


-100000E+01 -20U000E+Q1 .300000E+01 -400000E+01 -500000E+01 
-600000E+01 -700000€+01 


NItF IN OBS-CALC 


-00000 -.00000 .00000 .00000 -.00000 -~.00000 -.00000 
DERIVAITVES OF EACH FARAMETER 
-.00002 -.00005 .14980 -.00000 
-.00002 63116 -,00000 -.00001 
Uae s 2 -,00UU00 -.00000 .00530 
-.00000 -.00000 .00287 -.00u00 
-.00000 00374 -.00000 -.00000 
-00110 
DIFF IN OBS-CALC 
-.10000 -.00(00 .00000 .00000 -.00000 -.00000 -.00000 
DERIVATIVES OF FACKH PARAMETER 
-.00002 -.000uU5 .143535 -.00000 
-.00002 S032 -.00700 -.00001 
cULT00 -.00000 -.00000 (0515 
-.00006 -~,00000 .00279 -.00000 
-.900000 Pi (sito, -.00000 -.00000 
~O0107 
DItF IN OBS-CALC 
-00009 -.00u0 .00000 .00000 -00000 -.00000 -.00000 
DE“ITVATITVES OF EACH F4RAME] ZH 
-.00002 -,canus .14929 -.00000 
-.0G0002 CS 021 -,00000 -.00901 
-01099 -.00060 -.00000 .00515 
-.00000 -,000U0 .00279 -.00000 


f2 


- .00000 .00167 -.00000 -.00100 


-00107 
DIFF IN OBS-CALC 
~ 30006 -., 00000 .d0000 .00v00 -.00000 -.00000 -.00000 


rg a a yi 


THE BEST VALUE> Lr Fe ee Are 


4 
-3620037633E20S .4458529789E+02 .53619582/79E-01 





ExTIMATES OF THE EFRRIR JN GACH PARAMCIER Awe «STANDARD DEVIATION) 


©1248909550E-U0 .3815987E74E-01 .1455978757E-04 


NUMBER OF TIERA!ITUNS = $ [4—E SUM OF THE SuJARES OF THE EXRORS JS NOW = .282 


OBSERVED 4INUS CALCUI ATED VALUES OF HE Be> CIT Ant 


»31061E-07 =.83954--06 .241576-05 .10739E-05 -.11070F-05 -.26591F-05 -. 34959 
STP 


TIME. 2 MINUTES AND 15 seCUNDS 


SA 
WN 


APPENDIX VI 
DIMENSIONLESS TIME, +t, DEPENDENCE ON REAL TIME, 9 


From the definition 


a een 3 (1) 


Lr SS 


and the specific column parameters given in Appendix III, the 
numerical equivalence of equation (1) can be calculated. The calcu- 


lation is as follows: 


og 7 Osos | evaluated at 4th stage from ! 


initial conditions 


Whe 


Using the ue from Appendix VII 
Pi = B01 6) 5 = ,4965426 (3) 


From Appendix III 


1s 
bs F 3 _ 
Pp j = Ry j fosizy oy 3= .4144 (ft) (4) 


yeaa 
Pi 
and 


WW, ss = 4,4228 (moles) 


Therefore, 
t = 3.748 6 
or 


1 7+ unit = 0.267 (minutes) 


94 


APPENDIX VII 


INITIAL AND FINAL STEADY STATE VALUES OF COLUMN VARIABLES 
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LIQUIO MOLE FRACTIONS LISTED AS COMPONENTS PER PLATE 
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